# Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures

a M G de Lima, a R da Silva, D A Rade, Noureddine Bouhaddi

► To cite this version:

a M G de Lima, a R da Silva, D A Rade, Noureddine Bouhaddi. Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures. Engineering Structures, 2010, 32, pp.1479 - 1488. 10.1016/j.engstruct.2010.01.028. hal-01511298

HAL Id: hal-01511298

https://hal.science/hal-01511298

Submitted on 20 Apr 2017

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire **HAL**, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures

A.M.G. de Limaª, A.R. da Silvaª, D.A. Radeª, N. Bouhaddiᵇ

The effective design of viscoelastic dampers as applied to real-world complex engineering structures can be conveniently carried out by using modern numerical optimization and/or model updating techniques. However, the large number of exact evaluations of the cost functions, combined with the typically high dimensions of large finite element models of industrial structures incorporating viscoelastic materials, makes the numerical processes very costly, sometimes unfeasible. Those difficulties motivate the study reported herein, in which a general strategy to improve the standard condensation methods by taking into account a priori information of the modifications into the viscoelastic zones is introduced. The proposed method can be used with any condensation procedure, including direct reductions and component mode synthesis.

**Keywords:**

Viscoelastic damping

Robust condensation

Component mode synthesis

## 1. Introduction

The use of viscoelastic materials has been regarded as a convenient strategy in many types of industrial applications, where these materials can be applied either as discrete devices or surface treatments at a relatively low cost [1–3]. However, it presents some inherent drawbacks such as the influence of operational and environmental factors (frequency, temperature, pre-loads, etc.) [1]. Also, viscoelastic dampers (specially, surface treatments) are prone to induce considerable mass additions. A typical example is the application of viscoelastic constrained layer patches in compressor systems. In fact, having selected a viscoelastic material with optimum damping capability over the desired temperature and frequency range of compressor's operation, higher damping is achieved by increasing the viscoelastic treatment that results in an augmentation of the mass of the compressor [4]. This last feature leads to the necessity of performing optimization aimed at achieving the desired performance while complying with design and construction constraints. However, for finite element (FE) models of engineering structures incorporating viscoelastic materials composed by many thousands of degrees-of-freedom (DOFs) (for example an aerospace structure can be composed by more than $10^6$ DOFs) the time to compute the exact evaluations during the iterative processes, performed on the full FE matrices, can become prohibitive [5,6]. Those difficulties motivate the study reported herein

with the intent to propose a general strategy to reduce the size of the viscoelastic structures for mainly several reasons: the cost of computation when one wants to extract the eigensolutions or to predict the behavior of the full structure, the optimization or model updating procedures and reliability-based optimization which require fast iterative techniques.

The difficulty of any model reduction procedure of viscoelastic systems lies on the fact that the viscoelastic stiffness matrix depends on frequency and temperature. Moreover, supplementary difficulty is to complete the representation basis in order to reduce truncation effects. One strategy is to generate a set of Ritz vectors capable of representing with precision the structural behavior under a wide variety of structural modifications. For example, Balmès and Germès [6] studied the possibility of using a constant Ritz basis to create parametric families of reduced models. Bouazzouni et al. [7] developed an optimal method to construct additional vectors by using the dynamic behavior of the structure before modification. The disadvantage of such approaches lies in the fact that a basis of reduction composed by a great number of residual static vectors is obtained, augmenting the computation effort involved in the condensation.

The substructuring or also known as Component Mode Synthesis (CMS) method currently plays a considerable role in the analysis of the viscoelastic systems, where a geometrically complex structure can be subdivided into components called substructures or superelements, which are analyzed and condensed separately, while preserving the junction's DOFs between the substructures. The substructures are represented by their modes (in the sense of Ritz vectors), including the normal, the rigid body, the static, and

<footnote>ª Federal University of Uberlândia, School of Mechanical Engineering, Campus Santa Mônica, P.O.Box 593, CEP 38400-902, Uberlândia/MG, Brazil</footnote>

<footnote>ᵇ FEMTO-ST Institute, UMR 6174, R. Chaléat Applied Mechanics Laboratory, University of Franche-Comté, 24 Chemin de l'Epitaphe 25000, Besançon, France</footnote>
the interface modes. The structural modifications of a substructure make it possible to reanalyze the full system with a low cost of computation and make it possible to prepare each component independently of the others. This advantage enables the analysis of large-scale viscoelastic problems with the current data-processing resources.

Several CMS methods have been developed for forty years, and all of them make use of the vibration normal modes of the substructures. Depending on the boundary conditions applied to the substructure interfaces when these normal modes are obtained, the CMS methods can be classified into four groups: fixed interface methods [8–11]; free interface methods [10,11]; hybrid interface and loaded interface methods [11]. The variants of each group differ mainly in the choice of the supplementary Ritz vectors, and in the coupling procedure. In this paper one proposes an improvement of the classic Craig–Bampton Transformation (CBT) adapted for viscoelastic systems, by the introduction of residual static vectors associated to the external loads and the forces due to the viscoelastic damping.

In the remainder, after the various aspects related to the theoretical foundations, the description of a numerical application composed by a stiffened panel treated by a passive constraining damping layer demonstrates the effectiveness of the robust densification strategy of viscoelastic systems.

## 2. Finite element models incorporating viscoelastic damping

Consider the following finite element equations of motion of a viscoelastic structure in the frequency domain containing N DOFs [5]:

$$[\mathbf{K}_e + \mathbf{K}_v(\omega, T) - \omega^2 \mathbf{M}] \mathbf{Q}(\omega, T) = \mathbf{F}(\omega) \quad (1a)$$

$$F(\omega) = \mathbf{b} \mathbf{u}(\omega); \quad \mathbf{y}(\omega, T) = \mathbf{c} \mathbf{Q}(\omega, T) \quad (1b)$$

where $\mathbf{M} \in \mathbb{R}^{N \times N}$ is the mass (symmetric, positive-definite) matrix, and $\mathbf{K}_e \in \mathbb{R}^{N \times N}$ and $\mathbf{K}_v(\omega, T) \in \mathbb{R}^{N \times N}$ are the stiffness matrices (symmetric, nonnegative-definite) corresponding to the purely elastic and viscoelastic substructures, respectively. $\mathbf{Q}(\omega, T) \in \mathbb{R}^N$ and $\mathbf{F}(\omega) \in \mathbb{R}^N$ are, respectively, the vectors of displacements and external loads. $\mathbf{y}(\omega, T) \in \mathbb{R}^c$ is the vector of complex responses and $\mathbf{u}(\omega) \in \mathbb{R}^f$ is the reduced vector of external loads. Matrices $\mathbf{b} \in \mathbb{R}^{N \times f}$ and $\mathbf{c} \in \mathbb{R}^{c \times N}$ are Boolean matrices which enable to select, among the d.o.f.s, those in which the excitation forces are applied and responses are computed, respectively. In expressions (1), the other forms of damping and the dissipation of the original structure are neglected.

It is widely known that the dynamic behavior of viscoelastic materials depend on a number of factors, among which the most relevant are the excitation frequency and the temperature [1]. Various mathematical models have been developed to represent this behavior and have been shown to be particularly suitable to be used in combination with FE discretization [12–15]. In this paper, as the interest is confined to frequency-domain analyses, the so-named Complex Modulus is used in combination with the Frequency–Temperature Correspondence Principle and the Elastic–Viscoelastic Correspondence Principle. According to the Complex Modulus approach, the dynamic behavior of viscoelastic materials in the frequency domain can be represented by introducing frequency- and temperature-dependent complex material moduli as follows [1]:

$$G(\omega, T) = G'(\omega, T) [1 + i\eta_G(\omega, T)] \quad (2)$$

where $\omega$ and $T$ denote, respectively, the excitation frequency and the temperature of the viscoelastic material, $G'(\omega, T)$, $\eta_G(\omega, T) = G''(\omega, T)/G'(\omega, T)$ and $G''(\omega, T)$ designate, respectively, the so-named storage modulus, loss factor and loss modulus. It should be noted that definition (2) apply to both longitudinal and transverse moduli $E(\omega, T)$ and $G(\omega, T)$. However, it has been assumed

by most authors [5,6,12–15], that the Poisson ratio is independent from frequency and temperature in such a way that $E(\omega, T)$ and $G(\omega, T)$ are related to each other through the relation $G(\omega, T) = E(\omega, T)/[2(1 + \nu)]$. Such an assumption, which has been adopted in this study, has been argued by [16], who verified experimentally the variations of the Poisson ratio for PVC specimens.

The Frequency–Temperature Superposition Principle (FTSP), also known as Williams, Landell and Ferry (WLF) Principle [1], establishes an equivalence between the effects of the excitation frequency and temperature on the properties of a broad class of thermorheologically simple viscoelastic materials. This implies that the viscoelastic characteristics at different temperatures can be related to each other by changes (or shifts) in the actual values of the excitation frequency. This leads to the concepts of shift factor and reduced frequency, symbolically expressed as:

$$\begin{aligned} G(\omega, T) &= G(\omega_r, T_0) = G(\alpha_T \omega, T_0), \\ \eta(\omega, T) &= \eta(\omega_r, T_0) = \eta(\alpha_T \omega, T_0) \end{aligned} \quad (3)$$

where $T$ indicates an arbitrary value of the temperature, $T_0$ is a reference value of temperature, $\omega_r = \alpha_T(T)\omega$ is the reduced frequency, $\omega$ is the actual excitation frequency, and $\alpha_T(T)$ is the shift function. Fig. 1 illustrates the FTSP, showing that having the modulus and loss factor of an arbitrary viscoelastic material for different temperature values, $T_{-1}$, $T_0$, $T_1$, if horizontal shifts along the frequency axis are applied to each of these curves, all of them can be combined into a single one, named master curve. The horizontal shift is given by $\alpha_T$ and depends on the temperature.

Functions $G(\omega_r)$ and $\alpha_T(T)$ can be obtained from experimental tests for specific viscoelastic materials [1]. Drake and Soovere [17] suggest analytical expressions for the complex modulus and shift factor for various commercial viscoelastic materials. Eqs. (4) represent the complex modulus and shift factor as functions of temperature and reduced frequency in the intervals $210 \le T \le 360$ K and $1.0 \le \omega \le 1.0 \times 10^6$ Hz, for the 3M™ ISD112 viscoelastic material, as provided by those authors. The 3M™ ISD112 is a rubber-like viscoelastic polymer which is provided by the manufacturer in the form of adhesive tapes [18].

$$G(\omega_r) = 430700 + \frac{1200 \times 10^6}{1 + 3.241 \times (i\omega_r/1543000)^{-0.18} + (i\omega_r/1543000)^{-0.6847}} \quad (4a)$$

$$\begin{aligned} \log(\alpha_T) = & -3758.4 \times \left(\frac{1}{T} - \frac{1}{290}\right) + 225.06 \\ & \times \log\left(\frac{T}{290}\right) + 0.23273 \times (T - 290) \end{aligned} \quad (4b)$$

where $T$ is the temperature in Kelvin, $\omega_r$ is the reduced frequency in rad/s, and $G$ is the complex modulus in $N/m^2$.

Fig. 2 depicts the standardized curves representing the variations of the storage and loss moduli, the loss factor and the shift factor as functions of the reduced frequency and temperature, respectively, as obtained from Eqs. (4a) and (4b).

According to the Elastic–Viscoelastic Correspondence Principle [19] the derivation of the FE model accounting for the viscoelastic behavior can be carried out in two distinct phases: first, the element and global stiffness matrices are obtained by considering pure elastic behavior (hence, frequency- and temperature-independent material moduli), accounting for the strain state assumed; then, the material moduli are modified to account for the viscoelastic behavior (according to the model expressed by Eq. (2), in the present case). This implies that $G(\omega, T)$ can be factored-out of the matrix $\mathbf{K}_v(\omega, T) = G(\omega, T)\bar{\mathbf{K}}_v$ (where $\bar{\mathbf{K}}_v$ is a constant matrix), which can be combined with Eqs. (1) to produce the recep-tance or frequency response function (FRF) matrix:

$$\mathbf{H}(\omega, T) = \mathbf{c}\mathbf{Z}(\omega, T)^{-1}\mathbf{b} \quad (5)$$
<figure><img src="image_2.png" /><figcaption>Fig. 1. Illustration of the Frequency-Temperature Superposition Principle (FTSP).</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Master curves (—G; - - - η) for the 3M™ ISD112 viscoelastic material.</figcaption></figure>

where $Z(\omega, T) = K_e + G(\omega, T)\bar{K}_v - \omega^2 M$ is the so-named dynamic stiffness matrix.

The difficulty in predicting the vibration response in the time domain and performing eigenvalue analysis for the viscoelastic systems comes from the fact that the viscoelastic stiffness matrix depends on frequency and temperature. As a result, one has a nonlinear eigenvalue problem that must be solved iteratively [20]. Some procedures for dealing with this problem have been suggested, such as the recent contributions by Palmeri and Ricciardelli [21], in which the eigenvalue problem of viscoelastic systems has been derived in the time domain by the concept of modal relaxation functions. Yuan and Agrawal [22] and by Wagner and Adhikari [23] proposed an alternative state-space approach for the time-domain analysis of viscoelastic systems. Other alternatives have been proposed based on the adoption of

particular representations for the frequency-dependent behavior of the viscoelastic materials, such as the Fractional Derivative [12], Golla-Hughes-McTavish [13,14] and Anelastic Displacement Field [15] models, which enable to transform the equations of motion of the viscoelastic systems 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.

The interest herein is focused on frequency-domain responses, and with this aim, Eq. (5) can, in principle, be directly used for calculating the steady-state harmonic responses of the viscoelastic structures in the frequency domain. However, such a procedure can be unfeasible when large-scale finite element models are dealt with. In this case, the computations can be alleviated by using model reduction techniques.
### 3. Robust condensation of viscoelastic systems for design procedure

In the case of complex engineering structures of industrial interest, FE models are usually constituted by a large number of DOFs (hundreds of thousand or even millions). In such cases, it becomes practically impossible to compute the FRFs directly from Eq. (5), owing to the prohibitive computation times and storage memory required. This fact motivates the use of model reduction procedures, which aims at reducing the model dimensions (and the associated computational burden), 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 Eqs. (1), can be approached by projections on a reduced vector basis as follows:

$$Q(\omega, T) = T\hat{Q}(\omega, T) \qquad (6)$$

where $T \in C^{N \times NR}$ is the transformation matrix formed column-wise by a vector basis, $\hat{Q}(\omega, T) \in C^{NR}$ are generalized coordinates, and $NR \ll N$ is the number of reduced vectors in the basis. The generalized coordinates representing the contribution of each column of $T$ are chosen arbitrarily in which the reduce model provides a reasonable predictive capacity into a frequency bandwidth. Frequency band of interest is taken into account by computing a number of normal modes and retaining those below a certain frequency (1.5 times the last frequency of interest is typically).

By considering Eqs. (1) and (6), the transfer function (5) can be rewritten as:

$$\hat{H}(\omega, T) = c\hat{Z}(\omega, T)^{-1}\mathbf{b} \qquad (7)$$

where $\hat{Z}(\omega, T) = T^T K_e T + G(\omega, T) T^T \bar{K}_v T - \omega^2 T^T M T$. $\hat{Z}(\omega, T) \in R^{NR \times NR}$ is the reduced dynamic stiffness matrix that must be computed frequency by frequency and inverted in a direct way by using efficient numerical algorithms. For models containing viscous or structural damping, it is relatively common to use constant projection basis formed by the eigenvectors of the associated conservative structure, as the mass and stiffness matrices are invariant [1]. However, for viscoelastic systems, the selection of the reduction basis is more delicate as this condition does not hold. Owing to the dependence of the stiffness matrix with respect to frequency and temperature, the reduction basis should be able to represent the changes of the dynamic behavior as frequency and temperature vary. Three procedures have been adopted regarding the computation of the reduction basis: (i) one can simply neglect this dependence by considering the stiffness matrix as being constant [6]. In this case, the reduction basis is also constant; (ii) one can use a constant reduction basis obtained by the resolution of the nonlinear eigenvalue problem associated to a frequency- and temperature-dependent stiffness matrix [20]; (iii) one can use an iterative method for the re-actualization of the reduction basis according to frequency [6].

In this work, the strategy proposed consists in using a reduction basis formed by a constant modal basis (named herein nominal reduction basis), enriched by static residual vectors or equivalently static responses to account for the static effects of the modal truncation. These static responses are computed by using the tangent stiffness matrix representing the static behavior of the viscoelastic materials. As can be seen in Fig. 2(a), in the low frequency range, as the modulus and loss factor curves are prolonged by asymptotes, the extrapolation leads to the real values $G_0$ and $\eta_0 = 0$. On the other hand, for high frequencies, the extrapolation gives the complex values $G_\infty$ and $\eta_\infty$. The tangent stiffness matrix can thus be calculated as follows [6]:

$$K_0 = K_e + G_0 \bar{K}_v. \qquad (8)$$

The nominal basis can be obtained by the resolution of the eigenvalue problem:

$$\begin{gathered} (K_0 - \lambda_i M) \phi_i = 0 \quad i = 1, \dots, N \\ \phi_0 = [\phi_1 \ \phi_2 \ \dots \ \phi_{NR}], \quad \Lambda_0 = \text{diag}(\lambda_1, \dots, \lambda_{NR}). \end{gathered} \qquad (9)$$

This basis is enriched by introducing the residues formed by the static displacements associated to external forces,

$$R = K_0^{-1}\mathbf{b} \qquad (10)$$

that must be further completed by the residual vectors associated to the viscoelastic damping forces:

$$R_v^0 = K_0^{-1} \bar{K}_v \phi_0. \qquad (11)$$

As detailed in Reference [6], these residuals are interpreted as the columns of the flexibility matrix of the associated undamped system, associated to the coordinates of application of two types of forces to it: the external excitation forces and the damping forces. These latter can be better understood by examining (1), noting that the term involving the viscoelastic behavior can be moved to the right-hand side, where it plays the role of additional forces applied to the associated conservative structure. Thus, the enriched basis of reduction for the viscoelastic system is given as follows:

$$T_0 = [\phi_0 \ R \ R_v^0]. \qquad (12)$$

It is important to note that the time required to compute the exact matrix $K_0^{-1}$ for FE models composed by a large number of DOFs can become prohibitive. The numerical technique considered to be more accurate and time-effective consists in perform the Cholesky decomposition into the product of a lower triangular matrix and its conjugate transpose easier to be inverted [24].

Experience has demonstrated that the nominal basis (12) can be used to reduce the viscoelastic systems with reasonable accuracy, but is not capable of representing the modifications of the dynamic behavior provoked by the modifications which must be introduced into the model during iterative optimization or model updating processes. This means that the basis (12) should, in principle, be updated successively to guarantee a satisfactory accuracy of viscoelastic model reduction while the model evolves during iterations. However, this process would involve costly computations. To cope with this difficulty, the strategy suggested herein consists in performing a further enrichment of the nominal basis $T_0$ by a set of residual vectors calculated based on the knowledge of the parametric modifications of the viscoelastic zones.

#### 3.1. Robust static residual vectors

The robust condensation strategy lies in its use of a priori information of the nature and localization of the potential design modifications, while keeping the amplitude of the modification as a variable [25]. In such a way, for viscoelastic systems, the objective is to construct a set of robust frequency-independent Ritz vectors, to complete the standard reduction basis (12). By considering the Eqs. (1) and (5), for the modified structural configuration on the viscoelastic zones, the dynamic equilibrium equation in the frequency domain can be written as follows:

$$[\Delta Z_v(\omega, T) + Z(\omega, T)] Q(\omega, T) = F(\omega) \qquad (13)$$

where $\Delta Z_v(\omega, T) = \Delta K_v(\omega, T) - \omega^2 \Delta M_v$ is the variation of the dynamic stiffness matrix associated to the viscoelastic modifications. Eq. (13) can be interpreted as the dynamic equilibrium equation of the nominal viscoelastic model, subjected to the forces $F_\Delta(\omega, T)$ associated to the structural modifications:

$$Z(\omega, T)Q(\omega, T) = F(\omega) + F_\Delta(\omega, T) \qquad (14a)$$

where:

$$F_\Delta(\omega, T) = -\Delta Z_v(\omega, T)Q(\omega, T). \qquad (14b)$$

At this point, it is important to consider that, in the context of the present study, the modifications are introduced on a set of geometrical design parameters associated to the viscoelastic zones
that intervene in a rather complicated nonlinear fashion in those dynamic matrices. In an attempt to compute the modified matrices involved in the further condensation, it becomes interesting to perform a parameterization of the FE model, which is understood as a means of making them to appear explicitly in the finite element mass and stiffness matrices. This procedure enables to account for such modifications and/or uncertainties in the values of the design parameters in a straightforward way during iterative processes. Also, it facilitates, to a large extent, the evaluation of the sensitivities of the responses with respect to the design parameters [26–28]. In general, the parameters of mass and stiffness intervene nonlinearly into the FE matrices, and after manipulations the perturbed matrices having the design parameters factored-out can be expressed symbolically as follows:

$$ \Delta \mathbf{K}_v(\omega, T) = G(\omega, T) \bigcup_{i=1}^{\text{zones}} \bar{\mathbf{K}}_{vi}, \quad \Delta \mathbf{M}_v = \bigcup_{i=1}^{\text{zones}} \mathbf{M}_{vi} \qquad (15a) $$

$$ \bar{\mathbf{K}}_{vi} = \frac{\Delta p_i}{p_i} \sum_{\alpha} \alpha(p_i)^{\alpha} \bar{\mathbf{K}}_{vi}, \quad \mathbf{M}_{vi} = \frac{\Delta p_i}{p_i} \sum_{\beta} \beta(p_i)^{\beta} \mathbf{M}_{vi}. \quad (15b) $$

In the equations above, symbol $\bigcup$ indicates matrix assembling, subscript i indicate the modified zones (most frequently encompassing various viscoelastic elements) in which the parameter $p_i$ intervenes. $\bar{\mathbf{K}}_{vi}$ and $\mathbf{M}_{vi}$ designate the stiffness and mass matrices corresponding to those viscoelastic zones, respectively. These latter are decomposed according to Eq. (15b) by identifying the matrices $\alpha \bar{\mathbf{K}}_{vi}$ and $\beta \mathbf{M}_{vi}$ from which the $\alpha$- and $\beta$-order exponentials in the parameter $p_i$ can be factored-out.

As can be seen in Eq. (14b) the vector of forces associated to the structural modifications depends on the dynamic response of the modified structure, $\mathbf{Q}(\omega, T)$. Since this response is unknown a priori, those forces cannot be computed exactly. The essential concept of the robust condensation is expressed in the following [25]: (i) Eq. (14b) is used to generate a vector of forces which, even though it does not contain the exact forces associated to the modifications, will at least represent a subspace containing these vectors. This is accomplished by introducing the response of the nominal system into Eq. (14b); (ii) the resulting vector of forces is used to generate static responses, once again on the basis of the nominal model; (iii) the first two steps are repeated for each design parameter subjected to modifications. Many types of responses may be introduced into Eq. (14b), including the normal modes and the sensitivity vectors. If the vector introduced is composed by a truncated basis of normal modes, Eq. (14b) assumes the following form:

$$ \mathbf{F}_{\Delta}(\omega, T) \approx -\Delta \mathbf{Z}_v(\omega, T) \phi_0 \hat{\mathbf{Q}}(\omega, T). \qquad (16) $$

For example, for a parameter $p_i$ intervening in the viscoelastic mass and stiffness matrices, the basis of forces can be expressed under the following form:

$$ \mathbf{F}_{\Delta p_i} = \begin{bmatrix} \mathbf{F}_{\Delta p_i}^{\mathbf{M}_{vi}} & \mathbf{F}_{\Delta p_i}^{\bar{\mathbf{K}}_{vi}} \end{bmatrix} \qquad (17) $$

where $\mathbf{F}_{\Delta p_i}^{\mathbf{M}_{vi}} = \mathbf{M}_{vi} \phi_0 \Lambda_0$ and $\mathbf{F}_{\Delta p_i}^{\bar{\mathbf{K}}_{vi}} = G_0 \bar{\mathbf{K}}_{vi} \phi_0$ are the basis of forces associated to the viscoelastic mass and stiffness modifications. After obtaining the basis of forces, one can calculate a series of static responses of the nominal viscoelastic system based on the tangent stiffness matrix as follows:

$$ \mathbf{R}_{\Delta p_i} = \mathbf{K}_0^{-1} \mathbf{F}_{\Delta p_i}, \quad i = 1, 2, \dots, np. \qquad (18) $$

The final robust condensation basis taking into account a priori knowledge of the viscoelastic modifications can be expressed as follows:

$$ \mathbf{T} = [T_0 \quad \mathbf{R}_{\Delta}] \qquad (19) $$

where $\mathbf{R}_{\Delta} = [\mathbf{R}_{\Delta p_1} \quad \mathbf{R}_{\Delta p_2} \quad \dots \quad \mathbf{R}_{\Delta np}]$, and np indicates the number of design parameters to be modified. The residue matrix $\mathbf{R}_{\Delta}$ is not

<figure><img src="image_4.png" /><figcaption>Fig. 3. Block-diagram of a reanalysis process combining robust condensation.</figcaption></figure>

necessarily of maximum rank. Thus, with the aim of obtaining a limited number of independent residue vectors, it is appropriate to select the dominating directions of this basis, which can be done by performing the Singular Value Decomposition (SVD) of $\mathbf{R}_{\Delta}$ to identify its dominant singular values.

Fig. 3 illustrates a comparison between a cycle of optimization and/or model updating processes by using a standard reduction and the proposed robust condensation strategy. The robust strategy is used to approximate the behavior of the modified viscoelastic structures without the re-actualization of the nominal basis of reduction, leading to a drastic reduction of the time required for computing the responses. Moreover, this approach can be advantageously adapted to several other structural domains based on iterative processes: stochastic structural dynamics [28], nonlinear mechanics, and reliability-optimization-based design.

In the case of global modifications, for which the whole structure fully treated by viscoelastic materials is modified by using the same level of perturbation, the standard condensation procedure by using the reduction basis (12) is robust and do not need to be improved. In this case, there is only a frequency shift between the nominal model and the perturbed model, so that transformation (12) can represent correctly the perturbed viscoelastic model. But when the local modifications are introduced that is the case of structures partially treated by passive constraining layer damping, experience shows that the precision of this solution degenerates rapidly with the amplitude of the perturbation, demonstrating the interest in using the robust basis (19).

## 4. CMS combining robust condensation for viscoelastic systems

The objective of the substructuring methods is to calculate the dynamic behavior of a complex engineering structure, starting from the knowledge of the dynamic behavior of its substructures [8,9]. The principal differences between the substructuring methods are the choice of the basis of reduction for each substructuring, and the assembly technique between them. The principal sources of errors when implementing a substructuring method are: (i) the use of incomplete Ritz basis; (ii) errors in the discretization process of the geometry of an interface between two components; (iii) and the inequality of displacements and the balance of forces over an interface.

Fig. 4 illustrates two substructures finite element mesh whose DOFs can be partitioned in two subsets: an internal domain composed by all internal nodes belonging to the mesh (exception of the right edge), represented by the subscript i, and the nodes of the right edge defining the junction, represented by the subscript j.

The DOFs of the displacement vector $\mathbf{Q}(\omega, T)$ of a viscoelastic substructure SS2 are partitioned into two subsets: one related to the junction DOFS, $\mathbf{Q}_j(\omega, T)$, and the others associated to the interior DOFs, $\mathbf{Q}_i(\omega, T)$, according to the expression:

$$ \mathbf{Q}(\omega, T) = \begin{Bmatrix} \mathbf{Q}_j(\omega, T) \\ \mathbf{Q}_i(\omega, T) \end{Bmatrix}. \qquad (20) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Illustration of two substructures FE mesh with partitioned DOFs.</figcaption></figure>

The dynamic equilibrium equation of the substructure SS2 characterized by the tangent stiffness matrix (according to Eq. (8) in Section 3), can be written by taking into account the partition (20) as follows:

$$ \left( \begin{bmatrix} \mathbf{K}_{0jj} & \mathbf{K}_{0ji} \\ \mathbf{K}_{0ij} & \mathbf{K}_{0ii} \end{bmatrix} - \Lambda \begin{bmatrix} \mathbf{M}_{ij} & \mathbf{M}_{ji} \\ \mathbf{M}_{ij} & \mathbf{M}_{ii} \end{bmatrix} \right) \begin{Bmatrix} \mathbf{Q}_i \\ \mathbf{Q}_i \end{Bmatrix} = \begin{Bmatrix} \mathbf{F}_j \\ \mathbf{0} \end{Bmatrix} \quad (21) $$

where $Λ = diag(λ_1, ..., λ_{NS})$, NS is the number of DOFs of the substructure.

The assembly properties of the reduced mass and stiffness matrices depend on the type of selected modes and the assembly process. In the case of assembly level of reduced matrices of several viscoelastic substructures, in this paper, one uses the method suggested by Craig and Bampton [8], which imposes the continuity of displacements on the interface between the substructures, according to the following expression (for the system illustrated in Fig. 4):

$$ \mathbf{Q}_j^{(SS1)} = \mathbf{Q}_j^{(SS2)}(\omega, T). \quad (22) $$

In dynamic analysis of complex viscoelastic structures with a large number of DOFs, the Craig-Bampton can be used to drastically reduce the overall number of DOFs when combined with component mode synthesis technique. The condensation basis known as traditional Craig-Bampton Transformation (CBT) is composed by the vectors associated to the forced modes enriched by normal modes for a substructure with blocked junction. The constraint modes correspond to the static response of the substructure when a unit displacement is imposed on one selected DOF of the contour, and the others DOFs are blocked. For a viscoelastic substructure, they form the solution of the following equation:

$$ \begin{bmatrix} \mathbf{K}_{0jj} & \mathbf{K}_{0ji} \\ \mathbf{K}_{0ij} & \mathbf{K}_{0ii} \end{bmatrix} \begin{Bmatrix} \mathbf{I}_{jj} \\ \mathbf{G}_{ij} \end{Bmatrix} = \begin{Bmatrix} \mathbf{F}_j \\ \mathbf{0} \end{Bmatrix} \quad (23) $$

where $\mathbf{I}_{jj}$ represents a unit vector at DOF j, that is a vector with zero entries except at the DOF j that is equal to one. It represents a unit displacement imposed on the DOF j of the junction. Vector $\mathbf{G}_{ij}$ represents the static displacements due to zero forces at the interior DOFs. The resolution of Eq. (23) yields:

$$ \mathbf{G}_{ij} = (-\mathbf{K}_{0ii}^{-1} \mathbf{K}_{0ij}) \mathbf{I}_{jj}, \quad \mathbf{F}_j = (\mathbf{K}_{0jj} - \mathbf{K}_{0ji} \mathbf{K}_{0ii}^{-1} \mathbf{K}_{0ij}) \mathbf{I}_{jj}. \quad (24) $$

The normal modes correspond to the eigensolutions of the viscoelastic substructure. There are three different types of normal modes according to the boundary conditions imposed on the interface [9]. In this work, one uses the modes associated to the blocked interface, $\psi_i$, obtained by the resolution of the homogeneous viscoelastic problem associated to the internal DOF i, for which the DOF j is subjected to the boundary conditions. The vector representing the generalized coordinates are then composed by the physical coordinates $\mathbf{Q}_i$ and the modal coordinates $\mathbf{Q}_m$ associated to the normal modes with blocked interfaces $\psi$. The CBT transformation for a viscoelastic substructure is defined as follows:

$$ \mathbf{Q}_{SS2} = \begin{Bmatrix} \mathbf{Q}_i \\ \mathbf{Q}_i \end{Bmatrix} = \begin{bmatrix} \mathbf{I}_{ij} & \mathbf{0} \\ \mathbf{G}_{ij} & \psi \end{bmatrix} \begin{Bmatrix} \mathbf{Q}_i \\ \mathbf{Q}_m \end{Bmatrix} \quad (25a) $$

where:

$$ \mathbf{T}_{CBT} = \begin{bmatrix} \mathbf{I}_{jj} & \mathbf{0} \\ \mathbf{G}_{ij} & \psi \end{bmatrix}. \quad (25b) $$

The reduced problem is obtained by substituting Eq. (25b) into Eq. (1a) and premultiplying the resulted system equation by the transpose of the basis $\mathbf{T}_{CBT}$. However, for the same reason as explained in the previously section, for a viscoelastically damped substructure, the basis (25b) must be enriched by static residues associated to the external loads, and the forces associated to the viscoelastic modifications. The normal modes of the viscoelastic substructure chosen to construct the first order static residues are the modes associated to the blocked interface, and the residual vectors representing the static displacements related to the external forces, $\mathbf{R}$, according to Eq. (10), and the residues associated to the viscoelastic modifications with blocked junction, $\mathbf{R}_v^0$ and $\mathbf{R}_\Delta$, calculated according to Eqs. (11) and (18), respectively. Thus, the robust enriched CBT, $\mathbf{T}_{CBTE}$, for a viscoelastic substructure is:

$$ \mathbf{T}_{CBTE} = \begin{bmatrix} \mathbf{I}_{jj} & \mathbf{0} & \mathbf{0} & \mathbf{0} \\ \mathbf{G}_{ij} & \psi & \mathbf{R} & \mathbf{R}_v^0 \\ \mathbf{0} & \mathbf{0} & \mathbf{R}_v^0 & \mathbf{0} \end{bmatrix}^{\mathbf{T}_{CBTE}}. \quad (26) $$

The strategy for improving the robustness of the classical CBT basis (25b) with respect to structural viscoelastic modifications to obtain the basis (26), is similar to that described in Section 3.1, in which the optimal set of Ritz vectors $\mathbf{R}_\Delta$ are constructed based on the introduced modifications on the nominal viscoelastic model, to complete the standard enriched CBT basis $\mathbf{T}_{CBTE}^0$ for the viscoelastically damped substructure.

## 5. Numerical application

In this section, a numerical application is presented to illustrate the main features of the proposed robust condensation strategy of viscoelastic systems. Fig. 5 depicts the test structure composed by a freely suspended stiffened panel containing four stringers. The FE model without viscoelastic treatment is composed of 928 elements having a total number of 5940 DOFs. The viscoelastic treatment,

<figure><img src="image_6.png" /><figcaption>Fig. 5. FE model of the test structure (a); Substructure definition (b).</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** FRF amplitudes of the reference and reduced systems by using the nominal basis $T_{01}(a)$, $T_{02}(b)$ and $T_{03}(c)$.</figcaption></figure>

also indicated on Fig. 5(a), is composed of 10 viscoelastic patches, each one comprising 16 three-layer sandwich plate elements developed according to References [27,28]. The FE model of the partially viscoelastically treated panel contains 6840 DOFs. The geometric dimensions are: internal radius: 938 mm; length: 720 mm; arc length: 680 mm; thicknesses of the panel and the stringers: 1.5 mm and 0.75 mm, respectively; height of the stringers: 30 mm. The material properties for both panel and stringers are: Young modulus $E = 2.1 \times 10^{11}$ N/m$^2$; mass density $\rho = 7800$ kg/m$^3$; Poisson's ratio, $\nu = 0.3$. The material properties of the constraining layer are the same as those of the stiffened panel, and for the viscoelastic core, one uses the modulus function of the 3M™ ISD112 ($\rho = 950$ kg/m$^3$), as shown in Fig. 2(a). The nominal values of the thicknesses of the viscoelastic and constraining layers are 0.0254 mm and 0.5 mm, respectively.

## 5.1. Evaluation of the robust condensation basis

The interest here is to verify the robustness of the reduced viscoelastic model generated by the use of the proposed condensation method described in Section 3.1. The first test is intended to evaluate the nominal enriched basis of reduction by using the static residues associated to the external loads and the viscoelastic forces. The computations consisted in obtaining the driving point FRFs associated to point P indicated on Fig. 5(a), assuming a temperature of the viscoelastic material of 25 °C, and a frequency band of interest [135–210 Hz]. To verify the direct condensation, one considers the following nominal basis: $T_{01} = [\phi_0]$ (60 eigenvectors); $T_{02} = [\phi_0 \ R]$ (60 eigenvectors, plus one residual vector computed by Eq. (10)); $T_{03} = [\phi_0 \ R \ R_v^0]$ (60 eigenvectors, one residual vector computed by Eq. (10), 54 residual vectors computed according to Eq. (11) after SVD filtering). The residues $R_v^0$ were computed based on the largest singular values, for which the relation $\sigma_1/\sigma_i \le 1 \times 10^5$, for $i = 1$ to 60 holds.

Fig. 6(a), (b) and (c) show the FRFs computed by using the three nominal basis, as compared to the amplitudes of the response computed by using a reference basis formed by a far larger

number of eigenvectors (600) and residual vectors (600). The FRF amplitudes in [dB] have been computed by using a convenient reference factor through the relation Amplitude [dB] = $20 \times \log_{10}(|H(\omega, T)|/1e-6)$. It can be clearly seen that the accuracy is continuously improved upon successive enrichment of the reduction basis by the inclusion of residual vectors accounting for the static residues associated to the external loading and damping forces, to form the nominal basis $T_{02}$ and $T_{03}$, respectively. This is confirmed by the analysis of the relative errors between the approximations of the FRFs computed for the three basis of reduction, as shown in Fig. 7(a). Fig. 7(b) confirms that the use of first order residues associated with external forces and viscoelastic loads are sufficient to represent with accuracy the dynamic behavior of the viscoelastic damped system.

The interest now is to evaluate the robustness of the nominal basis further enriched to account for structural modifications introduced into the nominal model, according to Eq. (19). The modification considered consists in increasing the thickness of the constraining layer of the nominal system by 90%. Fig. 8 enables us to compare the amplitudes of the frequency response functions for the nominal and perturbed damped systems both computed using the same reference reduction basis, without further enrichment for structural modifications. It can be seen that the dynamic behavior of the perturbed viscoelastic damped system does not differ too much from that of the nominal system, although viscoelastic damping levels are strongly influenced.

In Fig. 9(a), the FRF of the perturbed reference system is compared to the counterpart computed by using the nominal basis $T_{03}$, containing 115 vectors. The observed differences lead to conclude that this basis is not capable of accurately representing the changes of the dynamic behavior induced by the structural modifications introduced. Fig. 9(b) enables us to compare the FRF of the perturbed reference system to the counterpart computed by using the basis $T_4 = [T_{03} \ R_\Lambda]$ containing 146 vectors (including 31 SVD-filtered residual vectors associated to the structural modifications). This time, one can observe a very satisfactory
<figure><img src="image_8.png" /><figcaption>Fig. 7. Comparison between the relative errors on the FRF amplitudes computed for the three reduction basis (a); relative errors computed for the basis $T_{03}$ (b).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. FRF amplitudes computed for the nominal and perturbed systems using the reference reduction basis.</figcaption></figure>

agreement between the amplitudes of the dynamic responses of both systems. This leads to conclude that the reduction basis $T_4$ is robust enough to represent the dynamic response of the perturbed viscoelastic structure in the frequency band of interest, yielding a

significantly more accurate prediction. Fig. 9(c) confirms that the use of residual vectors associated to the structural modifications is sufficient to represent with accuracy the dynamic behavior of the modified viscoelastic system.

## 5.2. Component mode synthesis

To verify the dynamic substructuring condensation strategy of viscoelastic structures, one uses the substructure definition represented in Fig. 5(b). The first substructure is composed by 551 shell finite elements, having 4350 DOFs (180 DOFs of junction), whose 160 finite elements are treated by 3M™ ISD112 viscoelastic material with the temperature of 25 °C, and the FE model of the second substructure is composed by a total of 377 finite elements, having 2520 DOFs, without viscoelastic treatment.

By using the enriched CBT basis (26), the nominal model of the first substructure is condensed to a substructure having 221 vectors described as follows:

$$T_{CBTE}^{0} = \begin{bmatrix} I_{180 \times 180} & 0_{180 \times 20} & 0_{180 \times 1} & 0_{180 \times 20} \\ G_{4170 \times 180} & \psi_{4170 \times 20} & R_{4170 \times 1} & R_{v_{4170 \times 20}}^{0} \end{bmatrix}$$

<figure><img src="image_10.png" /><figcaption>Fig. 9. FRF amplitudes computed for the nominal and perturbed systems using different reduction basis: (a) $T_{03}$; (b) $T_4$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. FRF amplitudes computed for the reference and reduced systems by using the enriched Craig–Bampton transformation.</figcaption></figure>

<table><caption>**Table 1**
Computing time for direct and CMS reduction methods.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Direct method</th><th colspan="2">CMS method</th></tr><tr><th><b><i>T</i><sub>01</sub> (s)</sub></th><th><b><i>T</i><sub>02</sub> (s)</sub></th><th><b><i>T</i><sub>03</sub> (s)</sub></th><th><b><i>T</i><sub>4</sub> (s)</sub></th><th><b><i>T</i><sub>CBTE</sub><sup>0</sup> (s)</sub></th><th><b><i>T</i><sub>CBTE</sub> (s)</sub></th></tr></thead><tbody><tr><td>Matrix assembly and reduction</td><td>18.3</td><td>18.4</td><td>19.9</td><td>3.3<sup>a</sup></td><td>14</td><td>3.8<sup>a</sup></td></tr><tr><td>Time to compute the reduced system (7)</td><td>2.8 s</td><td>2.9</td><td>3.7</td><td>5.8</td><td>17</td><td>19.9</td></tr><tr><td>Total computing time</td><td>21.1</td><td>21.3</td><td>23.6</td><td>9.1</td><td>31</td><td>23.7</td></tr></tbody></table>

$$ (221 \text{ vectors}) \qquad (27) $$

where $G_{4170×180}$ represents the constrained modes; $\psi_{4170×20}$ represents the normal modes associated to the blocked interface; $R_{4170×1}$ is the residues associated to the external loads; and $R_{4170×20}^0$ represents the residues associated to the viscoelastic forces.

In the same manner, the condensed model of the second substructure (purely elastic) is composed by 190 vectors, for which the classical CBT is expressed as:

$$ \mathbf{T}_{\text{CBT}} = \begin{bmatrix} \mathbf{I}_{180 \times 180} & \mathbf{0}_{180 \times 10} \\ \mathbf{G}_{2340 \times 180} & \psi_{2340 \times 10} \end{bmatrix} \quad (190 \text{ vectors}). \quad (28) $$

Fig. 10(a) shows the FRF amplitudes of the reference system and those obtained for the condensed model by using the enriched CBT. It is shown that the enrichment of the classical CBT basis of reduction associated to the external loads and the viscoelastic effects improves the results in the frequency band of analysis.

The modifications imposed on the nominal model consist in perturbing the viscoelastic zones of the damped substructure, by increasing the constraining layer thickness by 90%. The reduction of the modified system is made by using the nominal basis of reduction associated the nominal system (compound by 221 vectors), enriched by static residues evaluated on the modified matrices, to take into account the modifications on the viscoelastic parameters.

$$ \mathbf{T}_{\text{CBTE}} = [\mathbf{T}_{\text{CBTE}}^{0} \quad \mathbf{R}_{\Delta}] \quad (246 \text{ vectors}). \quad (29) $$

Fig. 10(b) enables us to compare the dynamic responses of the perturbed systems computed for the reference and reduced systems, respectively. Again, by this figure, one can conclude that the first order static correction associated to the viscoelastic modifications is robust enough to represent correctly the dynamic response of the perturbed system in the frequency band of interest.

To provide a sense of the additional computation effort associated with the enrichment process, Table 1 provides the comparison between the direct and CMS reduction methods. The most immediate use of the enrichment process is the increase of the total time required to compute the reduced system response (7). Also, the assembly of the matrices and their reduction take more time

of the total computing time for the direct reduction, which motivates the use of a *parameterization* of the FE model, leading to significant time savings in iterative processes. The results shown in Table 1 demonstrate that the CMS method is a little bit more costly than a direct method in term of computation but definitely more economic than an exact reanalysis when dealing with very large models having a large number of subassemblies. Moreover, one can take advantage of this enriched substructuring scheme when only one component incorporating surface viscoelastic treatment is subject to local modifications. In this case, there is no need for re-assembling or projecting on the reduction basis (27) the matrices of the others elastic components. Thus, the iterative processes would give better results in terms of computing time.

## 6. Concluding remarks

A robust condensation procedure combined with a substructuring technique, intended to be used for dealing with large-scale viscoelastically damped systems was suggested and evaluated. The original aspects of the procedure lies in the adaptation of the concept of robust condensation, initially developed for undamped structures, for systems containing viscoelastic materials, and the extension of the proposed methodology to the case of CMS method. This approach allows design parameters of a superelement to be modified, for example in the context of an optimization and/or model updating processes, without the necessity to perform a complete superelement analysis at each point in parameter space. Thus, the improved model reduction transformation can be prepared in advance for each substructure and then used directly during the iteration processes to avoid the exact recalculation of the modified zones.

For viscoelastic structures the proposed robust condensation methodology can be integrated in a variety of component mode synthesis techniques and we have illustrated its use in the context of the Craig–Bampton condensation procedure. Moreover, this approach is particularly effective when dealing with very large-scale FE viscoelastic models having a large number of substructures. At the present time, proposed design process is not constrained solely to the viscoelastic elements, and can already be used to approximate reanalysis of linear elastodynamic behavior to substructure technology.

<footnote>a Time to compute only the residues associated to the structural modifications and to perform the reduction.</footnote>
An academic example derived from the aeronautic engineer-
ing interest was used to illustrate the efficiency of the robust con-
densation procedure in representing the modified behavior of the
perturbed viscoelastic zones. The obtained results demonstrated
the effectiveness of the two robust condensation strategies mainly
in terms of the drastic reduction of the number of DOFs, which
demonstrates that the suggested technique is well adapted to be
applied to complex viscoelastic structures.

Acknowledgements

The authors wish to thank Brazilian Research Council – CNPq for the continued support to their research work, especially through research projects 480785/2008-2 (A.M.G. de Lima) and 310524/2006-7 (D.A. Rade); Minas Gerais State Agency FAPEMIG, CAPES Foundation, and French agency COFECUB for the support to the international cooperation project with the Applied Mechanics Laboratory R. Chaléat, FEMTO-ST.

References

[1] Nashif AD, Jones DIG, Henderson JP. Vibration damping. New York: John Wiley & Sons; 1985.

[2] Rao MD. Recent applications of viscoelastic damping for noise control in automobiles and commercial airplanes. In: Proceedings of the emerging trends in vibration and noise engineering symposium. 2001.

[3] Samali B, Kwok KCS. Use of viscoelastic dampers in reducing wind and earthquake induced motion of building structures. Eng Struct 1995;17(9): 639-54.

[4] de Lima AMG, Kotinda GI, Rade DA, Steffen V Jr, Baars E. Vibration attenuation in refrigerators using constrained viscoelastic layers: preliminary finite element modeling and experimental assessment. In: Proceedings of the 22th international conference on compressors and their systems. 2007.

[5] de Lima AMG, Rade DA, Lépore-Neto FP. An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring. J Mech Syst Signal Process 2009;23(4): 1272-81.

[6] Balmès E, Germès S. Tools for viscoelastic damping treatment design: Application to an automotive floor panel. In: Proceedings of the 28th international seminar on modal analysis. 2002.

[7] Bouazzouni A, Lallement G, Cogan S. Selecting a Ritz basis for reanalysis of frequency response functions of modified structures. J Sound Vib 1997;199(2): 309-22.

[8] Craig Jr RR, Bampton MCC. Coupling of substructures for dynamic analysis. AIAA J 1968;6:1313-9.

[9] Hurty WC. Dynamic analysis of structural systems using component modes. AIAA J 1965;3:678-85.

[10] MacNeal RH. A hybrid method of component mode synthesis. J Comput Struct 1971;4:581-601.

[11] Rubin S. Improved component mode representation for structural dynamic analysis. AIAA J 1975;13:995-1006.

[12] Bagley RL, Torvik PJ. A generalized derivative model for an elastomer damper. Shock Vib Bull 1979;49:135-43.

[13] Golla DF, Hughes PC. Dynamics of viscoelastic structures – A time domain finite element formulation. J Appl Mech 1985;52:897-906.

[14] McTavish D, Hughes PC. Modeling of linear viscoelastic space structures. J Vib Acoustics 1993;115:103-15.

[15] Lesieutre GA, Bianchini E. Time domain modeling of linear viscoelasticity using anelastic displacement fields. J Vib Acoustics 1995;117:424-30.

[16] Moreau A. Identification de propriétés viscoélastiques de matériaux polymères par mesures de réponses en fréquences de structures. Doctorate thesis. Rouen (France): INSA; 2007.

[17] Drake ML, Soovere J. A design guide for damping of aerospace structures, Proceedings of vibration damping workshop, AFWAL-TR-84-3064. In: Pub. by Air Force Wright Aeronautical Laboratories, Wright-Patterson Air Force Base, Ohio, 1984.

[18] http://www.3m.com/ (Accessed 03 August 2009).

[19] Christensen RM. Theory of viscoelasticity: An introduction. 2nd ed. New York: Academic Press Inc; 1982.

[20] Daya EM, Potier-Ferry M. A numerical method for nonlinear eigenvalue problems: Application to vibrations of viscoelastic structures. J Comput Struct 2001;79:533-41.

[21] Palmeri A, Ricciardelli F. Fatigue analyses of buildings with Viscoelastic dampers. J Wind Eng Indust Aerodyn 2006;94(5):377-95.

[22] Yuan L, Agrawal OP. A numerical scheme for dynamic systems containing fractional derivatives. J Vib Acoustics 2002;124(2):321-4.

[23] Wagner N, Adhikari S. Symmetric state-space method for a class of nonviscously damped systems. AIAA J 2003;41(5):951-6.

[24] Lay DC. Linear algebra and its applications. 3rd ed. Boston: Addison-Wesley; 2002.

[25] Masson G, Ait Brik B, Cogan S, Bouhaddi N. Component mode synthesis (CMS) based on an enriched Ritz approach for efficient structural optimization. J Sound Vib 2006;296(4):845-60.

[26] de Lima AMG, Faria AW, Rade DA. Sensitivity analysis of frequency response functions of composite sandwich plates containing viscoelastic layers. Compos Struct 2010;92(1):364-76.

[27] de Lima AMG, Stoppa MH, Rade DA, Steffen Jr V. Sensitivity analysis of viscoelastic structures. Shock Vib 2006;13(5):545-58.

[28] Guedri M, de Lima AMG, Bouhaddi N, Rade DA. Robust design of viscoelastic structures based on stochastic finite element models. Mech Syst Signal Process 2010;24(1):59-77.