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

HAL
open science

# A time-domain finite element model reduction method for viscoelastic linear and nonlinear systems

Antonio Marcos Gonçalves de Lima, Noureddine Bouhaddi, Domingos Alves

Rade, Marcelo Belonsi

► To cite this version:

Antonio Marcos Gonçalves de Lima, Noureddine Bouhaddi, Domingos Alves Rade, Marcelo Belonsi. A time-domain finite element model reduction method for viscoelastic linear and nonlinear systems. Latin American Journal of Solids and Structures, 2015, 12 (6), pp.1182-1201. hal-01512520

HAL Id: hal-01512520

https://hal.science/hal-01512520

Submitted on 23 Apr 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A time-domain finite element model reduction method for viscoelastic linear and nonlinear systems

Many authors have shown that the effective design of viscoelastic systems can be conveniently carried out by using modern mathematical models to represent the frequency- and temperature-dependent behavior of viscoelastic materials. However, in the quest for design procedures of real-word engineering structures, the large number of exact evaluations of the dynamic responses during iterative procedures, combined with the typically high dimensions of large finite element models, makes the numerical analysis very costly, sometimes unfeasible. It is especially true when the viscoelastic materials are used to reduce vibrations of nonlinear systems. As a matter of fact, which the resolution of the resulting nonlinear equations of motion with frequency- and temperature-dependent viscoelastic damping forces is an interesting, but hard-to-solve problem. Those difficulties motivate the present study, in which a time-domain condensation strategy of viscoelastic systems is addressed, where the viscoelastic behavior is modeled by using a four parameter fractional derivative model. After the discussion of various theoretical aspects, the exact and reduced time responses are calculated for a three-layer sandwich plate by considering nonlinear boundary conditions.

## Keywords

Passive vibration control; viscoelastic damping; fractional derivative model; condensation; nonlinear vibrations.

Antônio Marcos Gonçalves de Limaª

Noureddine Bouhaddiᵇ

Domingos Alves Radeᶜ

Marcelo Belonsiᵈ

ª,Federal University of Uberlândia, School of Mechanical Engineering, Campus Santa Mônica, Uberlândia-MG, Brazil

ᵇ,FEMTO-ST Institute, Chaléat Applied Mechanics Laboratory, Univeristy of Franche-Comté, Besançon, France

ᶜ,Instituto Tecnológico de Aeronáutica, ITA, São José dos Campos, SP, Brazil

Corresponding author:

ª,amglima@mecanica.ufu.br

ᵇ,noureddine.bouhaddi@univ-fcomte.fr

ᶜdomingos@ufu.br

ᵈ,marcelobelonsi@bol.com.br

## 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 (Nashif et al., 1985; Rao, 2001; Samali and Kwok, 1995). A typical example is the application of viscoelastic constrained layers in compressors systems (de Lima et al., 2007). However, the behavior of viscoelastic materials presents some inherent com-
plexities such as the influence of operational and environmental factors (frequency, temperature, pre-loads, etc.) (Nashif et al., 1985). Thus, several approaches have been developed for performing dynamic responses of engineering structures containing viscoelastic damping devices (Bagley and Torvik, 1979; Golla and Hughes, 1985; McTavish and Hughes, 1993; Lesieutre and Bianchini, 1995). However, many authors have found-out that the proposed mathematical models used for performing the frequency- and temperature-dependent behavior of viscoelastic materials, based on the concept of internal variables lead to global systems of equations of motion whose numbers of degrees of freedom (DOFs) largely exceed the order of the associated undamped structures (Balmès and Germès, 2002; de Lima et al., 2009; de Lima et al., 2010). As a result, if such response evaluations are made based on computations performed on the full finite element (FE) models of structures treated with viscoelastic materials composed by many thousands of DOFs, which is not rarely the case, the computational cost necessary to perform exact response evaluations during iterative processes such as structural optimization and uncertainty propagation, can become prohibitive, even unfeasible (de Lima et al., 2010; Zghal et al., 2014). Moreover, the increased dimension of the analytical models can preclude their use for active control and nonlinear vibration computations, in which time-domain analyses must be performed.

These drawbacks can be circumvented by combining directly the so-called model-reduction techniques in an attempt to reduce the order of the FE model while preserving its capability to represent the dynamic behavior of viscoelastic systems. However, it must be reminded that an inherent difficulty of any model reduction procedure applied to viscoelastic systems lies on the fact that the viscoelastic stiffness matrix depends on frequency and temperature. Another difficulty is to choose the projection bases in order to reduce truncation effects (de Lima et al., 2010).

One possible strategy is to generate a set of Ritz vectors capable of representing with satisfactory accuracy the structural behaviour under structural modifications. For example, Balmès and Germès (2002) studied the possibility of using a constant Ritz basis to create parametric families of reduced models developed for viscoelastic structures to be applied in the frequency-domain. Bouazzouni et al. (1997) developed an optimal method to construct additional vectors by using the dynamic behavior of undamped structures before modifications. De Lima et al. (2009) have proposed the use of a frequency- and temperature-independent reduction basis of viscoelastic systems based on the construction of residual static vectors taking into account the external loads and the viscoelastic damping forces.

The disadvantage of such approaches lies in the fact that a basis of reduction composed by a large number of residual static vectors is obtained, thus augmenting the computation effort involved in the condensation. Also, the proposed condensation strategies are restricted to the frequency-domain analysis, since the complex modulus approach is used to represent the viscoelastic dynamic features. More recently, Zghal et al. (2014) have proposed a model reduction method of nonlinear dynamic analysis of viscoelastic systems in time domain using the well-known Golla-Hughes-McTavish (GHM) model (Golla and Hughes, 1985; McTavish and Hughes, 1993). However, it leads to augmented global systems of equations of motion whose numbers of DOFs largely exceed the order of the associated undamped structures. This fact motivates the proposition of a time-domain condensation method based on the use of constant enriched reduction bases, which take into account a priori information of the viscoelastic damping forces and the local perturba-
tions to be applied in the transient analysis of viscoelastic systems. The perturbations can be applied related to either structural viscoelastic modifications or local nonlinearities.

As for the viscoelastic models, among the widely used mathematical representations accounting for the typical dependence of the viscoelastic properties with respect to frequency and temperature (Bagley and Torvik, 1979; Golla and Hughes, 1985; McTavish and Hughes, 1993; Lesieutre and Bianchini, 1995), in the present study, the proposed time-domain condensation approach is specially adapted for viscoelastic systems in which the viscoelastic behavior is modeled by using the four-parameter fractional derivative model (FDM) originally proposed by Bagley and Torvik (1979) and modified by Galucio et al. (2004). Although the association of this viscoelastic model in a FE discretization has been largely used in previous studies devoted to linear vibrating systems, the main contribution intended for the present study is the proposition of a straightforward time-domain condensation strategy to reduce linear and nonlinear viscoelastic systems, which requires specially adapted time-domain analytical and numerical resolution procedures.

In the remainder, after the presentation of various aspects related to the theoretical foundations, the description of numerical applications composed by a three-layer sandwich plate structure with embedded viscoelastic materials, subjected to linear and nonlinear structural modifications, demonstrates the effectiveness of the proposed condensation strategy.

## 2 REVIEW OF THE FE MODELING OF PLATES WITH VISCOELASTIC CONSTRAINING LAYERS

In this section, the model of a moderately thin three-layer sandwich plate FE, which can be frequently found, for example, in aerospace systems, is summarized based on the original developments made by Khatua and Cheung (1973) and implemented by de Lima et al. (2006). The inclusion of the frequency- and temperature- dependent behavior of the viscoelastic material is made by using the so-called *Elastic-Viscoelastic Correspondence Principle* (Nashif et al., 1985), according to which, the structural matrices are first generated for specific types of finite elements (rods, beams, plates, etc.) assuming that the longitudinal modulus and/or the shear modulus (according to the stress state assumed) are constant (independent on frequency and temperature). Then, after the finite element matrices are constructed, the frequency-temperature dependency of those *moduli* can be introduced according to the complex modulus approach combined with the *Frequency-Temperature Superposition Principle* (de Lima et al., 2009).

Figure 1 depicts a rectangular element formed by an elastic base-plate (1), a viscoelastic core (2) and an elastic constraining layer (3). This element contains four nodes and seven DOFs per node, representing the in-plane displacements in the middle plane of the base-plate in directions *x* and *y* (denoted by $u_1$ and $v_1$, respectively), the in-plane displacements of the middle plane of the constraining layer in directions *x* and *y* (denoted by $u_3$ and $v_3$, respectively), the transverse displacements, *w*, and the cross-section rotations about *x* and *y*, denoted by $\theta_x$ and $\theta_y$, respectively.

In the development of the theory, the following assumptions are adopted: (i) all the materials involved are homogeneous and isotropic and present linear behavior; (ii) normal strains in direction *z* are null for the three layers; (iii) the elastic layers (1) and (3) are modeled according to Kirchhoff's theory; (iv) for the viscoelastic core, Mindlin's theory is adopted (transverse shear is included); (v) the cross-section rotations, $\theta_x$ and $\theta_y$, are the same for the elastic layers; (vi) the
<figure><img src="image_3.png" /><figcaption>Figure 1: Schematic representation of three-layer sandwich plate element.</figcaption></figure>

transverse displacement, $w$, is the same for all the layers. These assumptions have been considered by many authors as being adequate for the modeling of thin panels, as is the case of the structures addressed herein (Austin, 1999). Moreover, previous studies carried-out by the authors demonstrated satisfactory correlation between model predictions and experimental results (de Lima et al., 2003).

The discretization of the displacement fields within the element is made by using bilinear interpolation functions for the displacements in the middle plane of the elastic layers in directions $x$ and $y$, and a cubic interpolation function for the transverse displacement, according to the relation, $\mathbf{u}(x, y, t) = \mathbf{N}(x, y) \mathbf{u}_{(e)}(t)$, where $\mathbf{N}(x, y)$ is the matrix formed by the shape functions, and $\mathbf{u}_{(e)}(t) = [u_1^i v_1^i u_3^i v_3^i w^i \theta_i^x \theta_i^y]^T$ with $i=1$ to 4 designates the vector containing the nodal variables as functions of time. The strain-displacement relations in vector form, $\boldsymbol{\varepsilon}(x, y, z, t) = \mathbf{B}(x, y, z) \mathbf{u}_{(e)}(t)$, are used and the resulting strains for the elastic layers, $\boldsymbol{\varepsilon}_k = [\varepsilon_x^{(k)} \varepsilon_y^{(k)} \gamma_{xy}^{(k)}]^T$ ($k=1,3$), and for the viscoelastic core, $\boldsymbol{\varepsilon}_2 = [\varepsilon_x^{(2)} \varepsilon_y^{(2)} \gamma_{xy}^{(2)} \gamma_{xz}^{(2)} \gamma_{yz}^{(2)}]^T$, are generated.

Thus, based on the stress-states assumed for each layer and the corresponding stress-strain relations, following standard analytical developments based on variational approaches, the strain and kinetic energies of the three-layer sandwich plate FE are formulated and the elementary mass and stiffnesses matrices are obtained as follows:

$$ \mathbf{M}^{(e)} = \sum_{k=1}^{3} \rho_k h_k \int_{x=0}^{b} \int_{y=0}^{a} \mathbf{N}^{T}(x, y) \mathbf{N}(x, y) \, dy \, dx \qquad (1.a) $$

$$ \mathbf{K}_{e}^{(e)} = \sum_{k=1,3} \int_{z=0}^{h_k} \int_{x=0}^{b} \int_{y=0}^{a} \mathbf{B}_{k}^{T}(x,y,z) \mathbf{C}_{k} \mathbf{B}_{k}(x,y,z) \, dy \, dx \, dz \qquad (1.b) $$

$$ \mathbf{K}_2^{(e*)} = \int_{z=0}^{h_2} \int_{x=0}^{b} \int_{y=0}^{a} \mathbf{B}_2^T(x, y, z) \mathbf{C}_2^* \mathbf{B}_2(x, y, z) \, dy \, dx \, dz \qquad (1.c) $$

where $\mathbf{C}_k$ ($k=1,3$) is the isotropic elastic material properties matrix, and $\mathbf{C}_2^* = \mathbf{C}_v^*$ contains the frequency- and temperature-dependent material properties for the viscoelastic core. Matrix $\mathbf{B}(x, y, z)$ is formed by differential operators appearing in the strain-displacement relations for each layer, as detailed in de Lima et al. (2006). Also, $b$ and $a$ designate, respectively, the dimension
sions of the plate element in directions x and y, and $h_k$ and $\rho_k$ are the thickness and the mass density of the k-th layer, respectively.

According to the theory of the sandwich plate finite element summarized above and assuming that the Poisson ratio is independent from frequency and temperature in such a way that the frequency- and temperature-dependent longitudinal and transverse moduli are related to each other through the relation, $E_v^* = 2G_v^*(1 + \nu)$, the design parameters of mass and stiffness of each layer can be factored-out of the elementary matrices by uncoupling membrane, bending and shear effects as:

$$ 
\begin{aligned}
\mathbf{M}_1^{(e)} &= \rho_1 h_1 {}_m\mathbf{M}_{1,1}^{(e)} + \rho_1 h_1 {}_b\mathbf{M}_{1,1}^{(e)} \\
\mathbf{M}_2^{(e)} &= \rho_2 h_2 {}_m\mathbf{M}_{2,0}^{(e)} + \rho_2 h_2 d_1 {}_m\mathbf{M}_{2,1}^{(e)} + \rho_2 h_2 d_1^2 {}_m\mathbf{M}_{2,2}^{(e)} + \rho_2 h_2 {}_b\mathbf{M}_{2,0}^{(e)} \\
\mathbf{M}_3^{(e)} &= \rho_3 h_3 {}_m\mathbf{M}_{3,1}^{(e)} + \rho_3 h_3 {}_b\mathbf{M}_{3,1}^{(e)}
\end{aligned} 
\quad (2.a) $$

$$ 
\begin{aligned}
\mathbf{K}_1^{(e)} &= E_1 h_1 {}_m\mathbf{K}_{1,1}^{(e)} + E_1 h_1^3 {}_b\mathbf{K}_{1,1}^{(e)} \\
\mathbf{K}_2^{(e*)} &= E_v^* h_2 \left( {}_m\mathbf{K}_{2,0}^{(e)} + d_1 {}_m\mathbf{K}_{2,1}^{(e)} + d_1^2 {}_m\mathbf{K}_{2,2}^{(e)} + {}_b\tilde{\mathbf{K}}_{2,0}^{(e)} + d_2 {}_b\mathbf{K}_{2,1}^{(e)} + d_2^2 {}_b\mathbf{K}_{2,2}^{(e)} \right) + \dots \\
\frac{E_v^*}{h_2} \left( {}_s\mathbf{K}_{2,0}^{(e)} + d_3 {}_s\mathbf{K}_{2,1}^{(e)} + d_3^2 d_2^2 {}_s\mathbf{K}_{2,2}^{(e)} \right) \\
\mathbf{K}_3^{(e)} &= E_3 h_3 {}_m\mathbf{K}_{3,1}^{(e)} + E_3 h_3^3 {}_b\mathbf{K}_{3,1}^{(e)}
\end{aligned} 
\quad (2.b) $$

where $d_3 = h_1 + 2h_2 + h_3$ and $E_v$ ($k=1,3$) represents the longitudinal moduli of the elastic layers. Also, subscripts m, b and s indicate, respectively, the membrane, bending and shear effects in the structural matrices.

It should be noted that, in Eqs. (2), the matrices appearing in the right-hand side are constant sparse matrices which are independent from the design parameters. It is clearly perceived how the use of such expressions facilitates a variety of finite element computations, since it enables to account for structural modifications and/or parametric 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.

From the elementary matrices computed for each element of the FE mesh, and neglecting other forms of damping, the elementary equations of motion are given as:

$$ \mathbf{M}^{(e)} \ddot{\mathbf{u}}_{(e)}(t) + [\mathbf{K}_{e}^{(e)} + \mathbf{K}_{v}^{(e*)}] \mathbf{u}_{(e)}(t) = \mathbf{f}^{(e)}(t) \quad (3) $$

where the stiffness matrices, $\mathbf{K}_e^{(e)} = \mathbf{K}_1^{(e)} + \mathbf{K}_3^{(e)}$ and $\mathbf{K}_v^{(e*)} = \mathbf{K}_2^{(e*)}$ give, respectively, the contributions of the purely elastic and viscoelastic parts of the damped structure, and $\mathbf{M}^{(e)} = \mathbf{M}_1^{(e)} + \mathbf{M}_2^{(e)} + \mathbf{M}_3^{(e)}$ includes the contribution of the base-plate, the viscoelastic core and the constraining layer to the inertia matrix. Finally, $\mathbf{u}_{(e)}(t)$ and $\mathbf{f}^{(e)}(t)$ are, respectively, the vectors of the generalized displacements and external loads.

Cleary, Eq. (3) must be used to construct the global equations of motion, accounting for the node connectivity of the discretization mesh, using standard FE assembling procedures.
## 2.1 Fractional derivative model incorporated into FE matrices

The dependency of the viscoelastic stiffness matrix on frequency and temperature is a consequence of the dependency of the material moduli with respect to these factors. A variety of mathematical models has been developed to represent the viscoelastic behavior and have been shown to be suitable to be used in combination with the FE discretization.

In this paper, as the interest is confined to a proposition of a time-domain reduction procedure of viscoelastic systems, the Fractional Derivative Model (FDM) initially proposed by Bagley and Torvik (1979) is used in combination with the Grünwald discretization technique (Galucio et al., 2004) to approximate the fractional operator appearing in the following one-dimensional stress state constitutive equation assumed for the viscoelastic material:

$$ \bar{\varepsilon}(t) + \tau^{\alpha} \frac{d^{\alpha} \bar{\varepsilon}(t)}{dt^{\alpha}} = \frac{E_{\infty} - E_0}{E_{\infty}} \varepsilon(t) \quad (4) $$

where $E_0$ and $E_{\infty}$ are, respectively, the relaxed or static modulus, and non-relaxed or high-frequency limit value of the modulus, $\tau$ is the relaxation time, and $\alpha$ represents the fractional order of the time derivative ($0 < \alpha < 1$).

One important aspect regarding the use of the following complex modulus function, $E_v^* = (E_0 + E_{\infty} (i\omega\tau)^{\alpha}) / (1 + (i\omega\tau)^{\alpha})$, obtained by applying the Fourier transform to Eq. (4) is the identification of the four parameters ($E_0, E_{\infty}, \tau, \alpha$) from experimental data-sheets provided by the manufacturers containing the material storage modulus and loss factor as functions of frequency and temperature. Thus, from the complex modulus definition, the determination of the values of the material parameters can be carried out by formulating an optimization problem in which the objective function represents the difference between the experimental data points and the corresponding model predictions in the frequency band of interest for a fixed temperature value (Galucio et al., 2004). As a result, the designer can obtain the storage modulus and the loss factor at any given temperature into a frequency band of interest, as illustrated in Fig. 2 for the particular material considered in the present study.

<figure><img src="image_4.png" /><figcaption>Figure 2: (a) Storage modulus and (b) loss factor for different temperatures for the 3M™ ISD112.</figcaption></figure>
Also, it must be emphasized that Eq. (4) was obtained by introducing the following internal variable as a strain function, $\bar{\epsilon}(t) = \epsilon(t) - \sigma(t)/E_{\infty}$, in the standard one-dimensional stress state constitutive equation for linear viscoelastic materials. As a result, Eq. (4) contains only one fractional derivative term, $d^{\alpha}\bar{\epsilon}(t)/dt^{\alpha}$, instead of two presented in the classical constitutive equation, and can be approximated by the following Grünwald definition, $d^{\alpha}\bar{\epsilon}(t)/dt^{\alpha} \approx \Delta t^{-\alpha} \sum_{j=0}^{n_p} A_{j+1}\bar{\epsilon}(t - j\Delta t)$, noting that $A_1 = 1$, to generate the following discretized form of the anelastic strain:

$$ \bar{\epsilon}(t) = (1-c) \frac{E_{\infty} - E_0}{E_{\infty}} \epsilon(t) - c \sum_{j=1}^{n_p} A_{j+1} \bar{\epsilon}(t - j\Delta t) \quad (5) $$

where $c = \tau^{\alpha}/(\tau^{\alpha} + \Delta t^{\alpha})$, $\Delta t = t/n$ is the time step, $n_p \le n$ is the number of discretization points, and $A_{j+1}$ represents the Grünwald coefficients given by the recurrence formula, $A_{j+1} = (j - \alpha - 1)A_j/j$.

At this point, the term $\mathbf{K}_v^{(e*)}\mathbf{u}_{(e)}(t)$ appearing in equations of motion (3) can be expressed as,

$$ \mathbf{K}_v^{(e*)}\mathbf{u}_{(e)}(t) = \int_z \int_y \int_x \mathbf{B}^T(x, y, z) \boldsymbol{\sigma}(x, y, z, t) \mathrm{d}x \mathrm{d}y \mathrm{d}z, $$

where for a general stress-strain relation, the vector $\boldsymbol{\sigma}(x, y, z, t)$ can be written in terms of the anelastic strain (5) taking into account the strain relation, $\bar{\epsilon}(t) = \epsilon(t) - \sigma(t)/E_{\infty}$, as follows:

$$ \boldsymbol{\sigma}(x, y, z, t) = E_0 \left[ \left( 1 + c \frac{E_{\infty} - E_0}{E_0} \right) \boldsymbol{\varepsilon}(x, y, z, t) - c \frac{E_{\infty}}{E_0} \sum_{j=1}^{n_p} A_{j+1} \bar{\boldsymbol{\varepsilon}}(x, y, z, t - j\Delta t) \right] \quad (6) $$

By combining Eq. (6) with the elastic and anelastic strain-displacement relations, $\boldsymbol{\varepsilon}(x, y, z, t) = \mathbf{B}(x, y, z)\mathbf{u}_{(e)}(t)$ and $\bar{\boldsymbol{\varepsilon}}(x, y, z, t) = \mathbf{B}(x, y, z)\bar{\mathbf{u}}_{(e)}(t)$, the resulting modified elementary equation of motion takes the following form:

$$ \mathbf{K}_v^{(e*)}\mathbf{u}_{(e)}(t) = \tilde{\mathbf{K}}_v^{(e)}\mathbf{u}_{(e)}(t) + \tilde{\mathbf{f}}_v^{(e)}(t - j\Delta t) \quad (7) $$

where

$$ \tilde{\mathbf{K}}_v^{(e)} = \left( 1 + c \frac{E_{\infty} - E_0}{E_0} \right) \mathbf{K}_v^{(e)} \quad \text{and} \quad \tilde{\mathbf{f}}_v^{(e)}(t - j\Delta t) = c \frac{E_{\infty}}{E_0} \mathbf{K}_v^{(e)} \sum_{j=1}^{n_p} A_{j+1} \bar{\mathbf{u}}_{(e)}(t - j\Delta t). $$

Upon introduction of Eq. (7) into Eq. (3) and after matrix assembling, the global equations of motion of the viscoelastic system containing N DOFs can be written as:

$$ \mathbf{M}\ddot{\mathbf{u}}(t) + (\mathbf{K}_e + \tilde{\mathbf{K}}_v)\mathbf{u}(t) = \mathbf{f}(t) - \tilde{\mathbf{f}}_v(t - j\Delta t) \quad (8) $$

where
$$ \boldsymbol{M} = \bigcup_{e=1}^{nel} \boldsymbol{M}^{(e)}, \quad \boldsymbol{K}_e = \bigcup_{e=1}^{nel} \boldsymbol{K}_e^{(e)} \quad \text{and} \quad \tilde{\boldsymbol{K}}^v = \bigcup_{e=1}^{nel} \tilde{\boldsymbol{K}}_v^{(e)} $$

are the global finite element mass and stiffness matrices, and symbol $\cup$ indicates matrix assembling. $\boldsymbol{u}(t)$ is the vector of global DOFs and $\boldsymbol{f}(t)$ is the vector of the generalized external loads. $\tilde{\boldsymbol{f}}_v(t - j\Delta t)$ is a load vector, which depends on the anelastic displacements, $\bar{\boldsymbol{u}}(t)$, and the Grünwald coefficients. These later account for the fading memory of viscoelastic materials.

It is interesting to highlight that Eq. (8) represents a time-domain model of structures containing viscoelastic materials, which can be interpreted as the result of incorporating modifications in the original elastic stiffness of the base-structure associated to viscoelastic materials, which is, in fact, independent of time, and by adding viscoelastic time-varying damping forces to the original external loads.

By comparing this modeling approach to other alternatives that have been proposed based on the adoption of particular representations for the frequency- and temperature-dependent behavior of the viscoelastic materials [5-8, 12], the proposed strategy does not require transformation of the equations of motion into a state-space form to generate frequency- and temperature-independent state matrices.

Based on Eq. (8), a time-domain condensation method based on the use of a frequency- and temperature-independent reduction bases is addressed next.

### 3 TIME-DOMAIN CONDENSATION OF VISCOELASTIC SYSTEMS FOR DESIGN PROCEDURE

In the case of complex 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 time-domain response directly from Eq. (8), owing to high computation times and storage memory required. This fact motivates the use of model reduction procedures, which aim at reducing the model dimensions (and the associated computational burden), while keeping a reasonable predictive capability of the numerical models. This can be done based on the assumption that the exact response, given by the resolution of Eq. (8), can be approached by a projection on a reduced vector basis as:

$$ \boldsymbol{u}(t) = \boldsymbol{T}_0 \hat{\boldsymbol{u}}(t) \qquad (9) $$

where $\boldsymbol{T}_0 \in C^{N \times NR}$ is the transformation matrix formed column-wise by vectors, $\hat{\boldsymbol{u}}(t) \in C^{NR}$ is the vector of generalized coordinates, and $NR \ll N$ is the number of vectors in the basis.

Once associated to the transformation expressed in Eq. (9), the equations of motion (8) can be rewritten as follows:

$$ \hat{\boldsymbol{M}} \ddot{\hat{\boldsymbol{u}}}(t) + (\hat{\boldsymbol{K}}_e + \hat{\boldsymbol{K}}_v) \hat{\boldsymbol{u}}(t) = \hat{\boldsymbol{f}}(t) - \hat{\boldsymbol{f}}_v(t - j\Delta t) \qquad (10) $$

where $\hat{\boldsymbol{K}}_e = \boldsymbol{T}_0^T \boldsymbol{K}_e \boldsymbol{T}_0$, $\hat{\boldsymbol{K}}_v = \boldsymbol{T}_0^T \tilde{\boldsymbol{K}}_v \boldsymbol{T}_0$, $\hat{\boldsymbol{M}} = \boldsymbol{T}_0^T \boldsymbol{M} \boldsymbol{T}_0$, $\hat{\boldsymbol{f}}_v(t - j\Delta t) = \boldsymbol{T}_0^T \tilde{\boldsymbol{f}}_v(t - j\Delta t)$ and $\hat{\boldsymbol{f}}(t) = \boldsymbol{T}_0^T \boldsymbol{f}(t)$ are the reduced matrices and vectors.
For models of structures containing viscous or structural damping, it is relatively common to use a constant projection basis formed by a subset of eigenvectors of the associated conservative structure, as the mass and stiffness matrices are invariant (Nashif et al., 1985). However, for systems containing viscoelastic materials, 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 (de Lima et al., 2010). This motivates the proposition of a special time-domain condensation method based on the use of a frequency- and temperature-independent enriched reduction basis which takes into account a priori information of the viscoelastic damping forces. These static responses are computed by using the tangent viscoelastic stiffness matrix computed in the low frequency range by considering the real value of the modulus, $E_0$, as follows:

$$ \mathbf{K}_0 = \mathbf{K}_e + E_0 \tilde{\mathbf{K}}_v \quad (11) $$

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

$$ \begin{gathered} (\mathbf{K}_0 - \lambda_i \mathbf{M}) \boldsymbol{\varphi}_i = 0 \quad i = 1, \dots, N \\ \boldsymbol{\varphi}_0 = [\boldsymbol{\varphi}_1 \ \boldsymbol{\varphi}_2 \ \dots \ \boldsymbol{\varphi}_{NR}], \quad \mathbf{A}_0 = \mathrm{diag}(\lambda_1, \dots, \lambda_{NR}) \end{gathered} \quad (12) $$

The basis, $\boldsymbol{\varphi}_0$, is further enriched by introducing the following residues formed by the static displacements associated, respectively, to viscoelastic damping forces and external excitations:

$$ \mathbf{R}_v^0 = \mathbf{K}_0^{-1} \tilde{\mathbf{K}}_v \boldsymbol{\varphi}_0 \quad (13.a) $$

$$ \mathbf{R} = \mathbf{K}_0^{-1} \mathbf{b} \quad (13.b) $$

where $\mathbf{b} \in \mathbb{R}^{N \times f}$ is a Boolean matrix which enables to select, among the DOFs, those in which the unit excitation forces are applied.

The residues (13.a) are interpreted as the columns of the flexibility matrix of the associated undamped system, related to the coordinates of application of the viscoelastic damping forces, which can be better understood by examining Eq. (8), 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:

$$ \mathbf{T}_0 = [\boldsymbol{\varphi}_0 \ \mathbf{R} \ \mathbf{R}_v^0] \quad (14) $$

## 3.1 Static residual vectors accounting for local nonlinearities

Local nonlinear behavior can be frequently found in a number of real-word engineering systems such as those having mechanical connections and joints. Thus, the objective of this section is to
extend the proposed reduction method to viscoelastically damped systems supported by lumped nonlinear springs. In this case, the inclusion of nonlinear effects into the global equations of motion (8) can be easily done by the concept of dyadic structural modifications (Maia and Montalvão e Silva, 1997), in which the loading system exerted by nonlinear springs on the viscoelastic system at the node of attachment, nd, can be expressed in the following vector form, by assuming cubic nonlinearity:

$$f_{nl}^{nd}(t) = K_{l}^{nd} u_{nd}(t) + K_{nl}^{nd} \tilde{u}_{nd}(t) \quad (15)$$

where $u_{nd}(t) = [u_1 \ v_1 \ u_3 \ v_3 \ w \ \theta_x \ \theta_y]^T$ and $\tilde{u}_{nd}(t) = [u_1^3 \ v_1^3 \ u_3^3 \ v_3^3 \ w^3 \ \theta_x^3 \ \theta_y^3]^T$ are the vectors containing the DOFs associated to the node of attachment of the linear and nonlinear contributions of the nonlinear springs in the global coordinate system.

Hence, the global system of equations of motion for the nonlinear viscoelastic system can be expressed under the following form:

$$M \ddot{u}(t) + (K_e + K_l + \tilde{K}_v) u(t) = f(t) - \tilde{f}_v(t - j\Delta t) - f_{\Delta}(t) \quad (16)$$

where $K_l = I_u^T K_l^{nd} I_u$ and $K_{nl} = I_{\tilde{u}}^T K_{nl}^{nd} I_{\tilde{u}}$ are, respectively, the diagonal stiffnesses matrices corresponding to the linear and nonlinear contributions, $I_u$ and $I_{\tilde{u}}$ designate the u, $\tilde{u}$-th columns of the identity matrices of order N, and $f_{\Delta}(t) = K_{nl} \tilde{u}(t)$.

For small nonlinearities, Eq. (16) shows that the nonlinear behavior can be understood as a perturbation indicated by, $f_{\Delta}(t)$, of the linear behavior, and can be viewed as a modification introduced on the linear viscoelastic system. As a result, the reduction method presented in the previous section can be extended to the case of nonlinear viscoelastic systems by a further enrichment of the reduction basis including normal modes related to the linear conservative associated system perturbed by the nonlinear springs. These modes can be obtained by the resolution of the eigenvalue problem:

$$(\tilde{K}_0 - \tilde{\lambda}_i M) \tilde{\phi}_i^0 = 0 \quad i = 1, \dots, N \quad (17)$$

where $\tilde{K}_0 = K_0 + K_l$.

Based on the fact that the vibration modes of the perturbed system does not differ appreciably from that of the non-perturbed nominal system, it is possible to assume that, $\tilde{\phi}_i = \tilde{\phi}_i^0 + \Delta\tilde{\phi}_i$, where, $\tilde{\phi}_i^0$, designates the modes of the nominal conservative associated system obtained by the resolution of Eq. (17), and $\Delta\tilde{\phi}_i$ is the perturbation due to the local nonlinearities which can be interpreted as an external load applied on the linear model, $\tilde{\phi}_i^0$. Thus, a set of static responses can be generated:

$$\Delta\tilde{\phi} \approx \tilde{K}_0^{-1} \tilde{b}_{\Delta} \quad (18)$$

where $\tilde{b}_{\Delta} \in R^{N \times m}$ is also a Boolean matrix which enables to select the nonlinear DOFs in which the unit forces are applied and static responses are computed, and m is the number of nonlinear DOFs.
Hence, the final condensation basis taking into account a priori knowledge of the perturbations can be expressed as follows:

$$ T = [T_0 \ \Delta\tilde{\phi}_1 \ \dots \ \Delta\tilde{\phi}_m] \qquad (19) $$

The residues, $\Delta\tilde{\phi}$, are not necessarily of maximum rank. Thus, with the aim of obtaining a limited number of independent residual vectors, it is appropriate to select the dominating directions of this basis, which can be done by performing the Singular Value Decomposition (SVD) of $T$ to identify its dominant singular values.

By considering the basis (19), the equations of motion (16) can be rewritten as:

$$ \hat{M} \ddot{u}(t) + \hat{K}_0 u(t) = \hat{f}(t) - \hat{f}_v(t - j\Delta t) - \hat{f}_\Delta(t) \qquad (20) $$

where $\hat{M} = T^T MT$, $\hat{K}_0 = T^T \tilde{K}_0 T$, $\hat{f}(t) = T^T f(t)$, $\hat{f}_v(t - j\Delta t) = T^T \tilde{f}_v(t - j\Delta t)$ and $\hat{f}_\Delta(t) = T^T f_\Delta(t)$ are the reduced matrices and vectors.

As will be demonstrated later on, the reduction basis given in Eq. (19) is robust enough to provide accurate predictions of the reduced model under the influence of the viscoelastic treatment, external excitations and nonlinear supports.

## 4 ITERATIVE RESOLUTION PROCEDURE

Based on the modeling features of the viscoelastic behavior and the presence of local nonlinearities, the resulting nonlinear differential equation of motion (20) were numerically solved using a variant of the Newmark integration scheme (Bathe, 1996) performed in the MATLAB environment.

Figure 3 shows the main steps of the numerical procedure which can be summarized as follows: ($\dot{i}$) at the beginning of the process, $t_0$, after defining the operating temperature of the viscoelastic material, $T_v$, the Newmark integration parameters, $(\beta, \gamma)$, and the variables $(E_\infty, E_0, \tau, \alpha)$ for FDM model obtained from curve-fitting of experimental data (Galucio et al., 2004), it is possible to compute the real and anelastic displacements, velocities and accelerations fields, $u(t_0), \bar{u}(t_0), \dot{u}(t_0)$ and $\ddot{u}(t_0)$, respectively, and the full mass, $M$, and stiffnesses matrices $(K_e, \tilde{K}_v, K_l, \tilde{K}_{nl})$ taking into account the mechanical boundary conditions; ($\dot{i}\dot{i}$) next, an eigenvalue problem is performed according to Eq. (17) in order to extract the eigenvalue modes, $\tilde{\phi}^0$, of the linear conservative associated system to generate the first basis of reduction to be further enriched by the residual static vectors associated to the external excitations, $R = \tilde{K}_0^{-1}b$, the viscoelastic damping forces, $R_v^0 = \tilde{K}_0^{-1} \tilde{K}_v \tilde{\phi}_0$, the modifications introduced on the viscoelastic zones, $\Delta\tilde{\phi} \approx \tilde{K}_0^{-1} \tilde{b}_\Delta$, or by the local nonlinearities, $\Delta\tilde{\phi} \approx \tilde{K}_0^{-1} \tilde{b}$. With the enriched basis of reduction it is possible to reduce the system in order to generate the reduced matrices, $\hat{M}$ and $\hat{K}_0$ and vectors, $\hat{u}, \hat{\bar{u}}, \hat{f}, \hat{f}_\Delta$ and $\hat{\tilde{f}}$; ($\dot{i}\dot{i}\dot{i}$) next, an iteration is initiated and the loading vectors due to the external loads, $\hat{f}$, the nonlinear perturbation, $\hat{f}_\Delta$, and the viscoelastic damping forces, $\hat{f}_v$, are evaluated at each iteration process. Also, a new set of reduced displacement fields, $\hat{u}_{i+1}$, are generated through the residues, $R_{i+1}$, performed according to the Newmark integration scheme, and taking into ac-
count the reduced matrices, $\hat{M}$ and $\hat{\mathbf{K}}_0$; (iv) in the next step, a new set of reduced velocities, $\dot{\hat{u}}_{i+1}$, and accelerations, $\ddot{\hat{u}}_{i+1}$, and the reduced anelastic displacements, $\hat{\mathbf{u}}_{i+1}$, fields, are generated. The iterative process is stopped when the maximum number of iterations is reached with a tolerance criterion supported by the Newmark integration scheme. Also, the iterative process is the same for linear viscoelastic systems, where $\hat{\mathbf{K}}_0 = \hat{\mathbf{K}}_e + E_0 \hat{\mathbf{K}}_v$ and $\hat{f}_{\Delta}(t) = 0$ are retained in Eq. (20).

<figure><img src="image_5.png" /><figcaption>Figure 3: Flowchart of the numerical procedure used to compute the system response.</figcaption></figure>

## 5 APPLICATIONS TO PLATE TREATED WITH CONSTRAINED VISCOELASTIC LAYER

In this first application, numerical tests were performed using the representative system shown in Fig. 4 composed by a clamped-clamped-free-free rectangular plate made of aluminum, fully treated with a layer of 3M ISD112™ (2013) viscoelastic material constrained between the base-plate and an outer aluminum sheet. The FE model is composed of 140 elements and the viscoelastic treatment is modeled according to the three-layer sandwich plate element developed according to Section 2. The geometric dimensions in x and y directions are depicted in the same figure, and the thicknesses of the base-plate, constraining layer, and viscoelastic core are, respectively, 1.0 mm, 0.5 mm and 0.254 mm. The material properties for the base-plate and constraining layer are: Young modulus $E = 70 \times 10^9$ N/m$^2$; mass density $\rho = 2700$ kg/m$^3$; Poisson ratio, $\nu = 0.3$. For the 3M™ ISD112, $\rho = 950$ kg/m$^3$, $\nu = 0.49$, and the parameters of the FDM model can be obtained by performing a curve-fitting procedure taking into account the viscoelastic material data provided in Fig. 2 in terms of storage modulus and loss factor, for a given temperature of the viscoelastic material.
The plate is subjected to a harmonic force of the form, $F(t) = F_0 \sin(2\pi f_0 t)$, where $F_0 = 10$ N is the amplitude of the excitation and $f_0$ is the excitation frequency, and to a transverse triangular impulse loading applied to point I indicated in Fig. 4. The numerical tests were performed in the time interval, $t \in [t_0 = 0$ s, $t_f = 0.35$ s], for an arbitrarily chosen constant time step, $\Delta t = 0.01$ ms, and the computations consisted in obtaining the normalized transverse displacements, $10^3 \times u(t)/L$, associated to point I.

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Sketch of the sandwich plate with the mechanical boundary conditions.</figcaption></figure>

First, the responses are computed using the full FE model without any condensation. Figures 5 represents the amplitudes of the normalized transverse displacements of the plate with and without viscoelastic treatment due to the applied impulse and harmonic loadings by assuming a temperature of the viscoelastic material of 26° C with the following FDM parameters, $E_0 = 1.2817$ MPa, $E_∞ = 454.52$ MPa, $τ = 5.94×10^{-7}$ s and $α = 0.6744$. For the harmonic loading, the excitation frequency ($f_0 = 76.47$ Hz) is close to the first flexural mode of the sandwich plate in which the viscoelastic stiffness matrix was computed in the low frequency range by considering the real value of the modulus, $E_0$. Also, it was considered a proportional damping computed by the following relation, $D = 0.005 \times K_e$. As expected, as the time increases, the flexural vibrations become larger showing the resonance characteristics of the plate for this loading condition. By comparing the time-domain responses of the sandwich plate for both input signals, it can be clearly perceived that a significantly reduction in the amplitudes of the transverse displacements of the plate is achieved, even when the excitation frequency is close to the first resonance.

The interest now is to evaluate the enriched basis of reduction by using the static residues associated to the external excitations and viscoelastic damping forces. To verify the direct condensation procedure, one considers the following nominal basis: $T_1 = [\varphi_0 \ R]$ (19 eigenvectors, plus one residual vector computed by the Eq. (13.b)); $T_2 = [\varphi_0 \ R \ R_v^0]$ (19 eigenvectors, one residual computed by the Eq. (13.b), 14 residual vectors computed according to the Eq. (13.a) after SVD filtering). The residues $R_v^0$ were computed based on the largest singular values, for which the relation, $\sigma_{max}/\sigma_i \le 1×10^6$, for $i = 1$ to 14 holds. The interest in examining these situations is to quantify the improvement entailed by the residual vectors associated to the viscoelastic
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Transverse displacements of the plate with and without surface viscoelastic treatment due to the applied harmonic (left) and impulse (right) loadings.</figcaption></figure>

damping forces and to demonstrate the capability of the time-domain reduction modeling procedure to accommodate such design variants.

Figure 6 shows the normalized time responses obtained by using the two nominal bases, as compared to the reference displacements computed by using the exact system. It can be clearly seen that the accuracy is improved upon enrichment of the reduction basis by the inclusion of residual vectors accounting for the static residues associated to the viscoelastic damping forces for both input signals, demonstrating that the use of these residues are sufficient to represent with accuracy the time-domain behavior of viscoelastically damped structures.

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Transverse displacements of the reference and reduced systems by using the nominal basis $T_1$ and $T_2$.</figcaption></figure>
As a complementary demonstration of the utility of the proposed condensation strategy in the analysis of viscoelastic systems, Fig. 7 shows the amplitudes of the frequency response functions (FRFs) of the sandwich plate computed by using the two nominal bases for the impulse loading condition, as compared to the amplitudes of the FRF computed for the exact system. Again, Fig. 7(b) confirms that the use of first order residues associated with viscoelastic damping forces are sufficient to represent with satisfactory accuracy the dynamic behavior of the viscoelastic system. Moreover, these results provide a sense of the efficiency of the surface viscoelastic treatment in mitigating the transverse vibrations of the plate.

<figure><img src="image_9.png" /><figcaption>Figure 7: Amplitudes of the FRFs of the reference and reduced systems by using the nominal basis $T_1$ and $T_2$.</figcaption></figure>

## 5.1 Structural modifications

In this second application the utility of the time-domain reduction strategy in the analysis of modified viscoelastic systems is shown. The interest is to evaluate the robustness of the nominal basis further enriched to account for small perturbations introduced into the nominal model, according to Eq. (19). The modifications considered consist in increasing the thicknesses of the viscoelastic and constraining layers of the nominal system in 80%, as indicated in Fig. 4, and the exact time-domain responses and the amplitudes of the FRFs of the perturbed system was computed using the modified FE model, as shown in Figure 8. It can be noted that the dynamic behavior of the perturbed system does not differ appreciably from that of the nominal system.

<figure><img src="image_10.png" /><figcaption>Figure 8: Time-domain responses and FRFs amplitudes for the nominal and perturbed exact systems.</figcaption></figure>
In Fig. 9, the responses of the perturbed exact system in the time and frequency domains are compared to the counterparts computed by using the nominal reduction basis $T_2$ containing 34 vectors, without further enrichment accounting for structural modifications. The observed differences for both domains lead to conclude that this basis is not capable of representing accurately enough the changes of the dynamic behavior induced by the structural viscoelastic modifications, since the basis $T_2$ represents with accuracy the behavior of the nominal system, as shown in Fig. 6. Thus, in a number of real-word engineering applications of passive constraining damping layers, this reduction may be found to be insufficient to provide the necessary representation of the dynamic behavior of viscoelastic materials associated to structural variations. In the same figure, it is possible to compare the dynamic responses of the perturbed exact system to the counterparts computed by using the basis $T_3 = [T_2 \ T_3]$ containing 41 vectors (including 7 SVD-filtered residual vectors) associated to the structural modifications. This time, one can observe a very satisfactory agreement between the amplitudes of the dynamic responses of both models. This leads to conclude that the reduction basis $T_3$ is robust enough to represent the response of the perturbed viscoelastic system.

<figure><img src="image_11.png" /><figcaption>**Figure 9:** Time-domain responses and FRFs amplitudes for the exact and reduced systems by using the basis $T_2$ and $T_3$.</figcaption></figure>
## 5.2 Modifications due to local nonlinearities

In this section, it is considered the same system shown in Fig. 4 subjected to non-linear boundary conditions through the application of transverse non-linear springs on the right edge of the plate located at $x = 0.30$ m. The non-linear behavior of the springs are computed by the Eq. (15), where the linear, $K_l = k_l \times \bar{K}_l$, and non-linear, $K_{nl} = k_{nl} \times \bar{K}_{nl}$, stiffnesses are evaluated by assuming the values $k_l = 100$ N/m and $k_{nl} = 1 \times 10^9$ N/m$^3$, where $\bar{K}_l$ and $\bar{K}_{nl}$ are diagonal identity matrices. For the harmonic loading condition, $F_0 = 30$ N is the amplitude of the excitation and the excitation frequency, $f_0 = 24$ Hz, is close to the first flexural mode of the conservative associated system in which the stiffness is computed by the relation, $\tilde{K}_0 = K_e + E_0 \tilde{K}_v + K_l$.

The interest is to evaluate the enriched basis of reduction by using the static residues associated to the non-linear perturbations by considering the following basis: $T_2 = [\varphi_0 \mathbf{R} \mathbf{R}_v^0]$ (19 eigenvectors, one residual vector computed by the Eq. (13.b), 14 residues computed according to the Eq. (13.b) after SVD filtering); $T_3 = [T_2 \Delta\tilde{\varphi}]$ (the basis $T_2$ and 11 vectors computed by the Eq. (18)).

Figure 10 shows that a reasonable agreement between the dynamic responses of both models is obtained by using the basis, $T_3$. The agreement is less satisfactory for the reduction basis, $T_2$, leading to conclude that the static residues associated to the local nonlinearities is robust enough to represent the response of the non-linear viscoelastic system.

<figure><img src="image_12.png" /><figcaption>Figure 10: Time-domain responses of the reference and reduced systems by using the basis $T_2$ and $T_3$ for 30 N.</figcaption></figure>

Clearly, it should be reminded that this kind of reasoning is strictly valid for small non-linear damped behaviors as, depending on the amplitude of the excitation, the use of a constant enriched basis can only provide rough approximations of exact non-linear behaviors, as shown in Fig. 11. It can be noted the degree of influence of the non-linear behavior on the prediction: the larger the non-linear behavior with respect to a given excitation, the larger the influence of it on the accuracy of the reduced model. In this case, the results could possibly be improved by updating the reduction basis, $T_3$, at each iteration cycle to account for the non-linear variations. In general, it results in a prohibitive computational time, since a set of updated eigenproblems must be solved. Another strategy is the so-called Combined Approximations Approach (CA) (Kirsch and Bogomolni, 2007) developed originally for linear and non-linear undamped reanalysis. These aspects should be considered in future works.
<figure><img src="image_13.png" /><figcaption>**Figure 11:** Time-domain responses of the reference and reduced systems by using the basis $T_2$ and $T_3$ for 50 N and 100 N.</figcaption></figure>

# 6 CONCLUDING REMARKS

A robust time-domain condensation procedure intended to be used for dealing with linear and non-linear 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 linear undamped structures in the frequency-domain, for systems containing viscoelastic materials in the time-domain, and the extension of the proposed methodology to non-linear structures incorporating viscoelastic materials with the aim of vibration attenuation. This approach allows design parameters of a superelement to be modified, for example in the context of linear and non-linear reanalysis such as an optimization and/or model updating processes and structural damage analysis, 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 and then used directly during the iteration processes to avoid the exact recalculation of the modified zones.

At the present time, the proposed time-condensation process is not constrained solely to the dynamic applications discussed herein, and can already be used to approximate reanalysis of linear and non-linear behaviors to viscoelastically damped structures in control technology and optimization and structural damage analyses where the iterative solutions consist of repeated analyses followed by redesign steps. However, it must be reminded that depending on the structural local perturbations introduced on the viscoelastic zones and the non-linear behavior to be considered, the proposed time-condensation are no longer valid and the reduction basis must be adapted to account for the local perturbations during the iteration cycles.

Finally, it is important to pointed-out that the efficiency of the proposed time-condensation compared with the complete analysis of the modified viscoelastic systems can be measured by the accuracy the dynamic responses, as demonstrated by the obtained results, and computational effort. It was found that the CPU effort for the dynamic reanalysis of the modified linear viscoelastic systems is reduced by more than 64%, compared with the complete analyses of the modified design. For the viscoelastic systems incorporating local non-linearities the total computational effort can be reduced by more than 96% for both cases investigated herein.
## Acknowledgements

The authors are grateful to the Minas Gerais State Agency FAPEMIG for the financial support to their research activities and the Brazilian Research Council – CNPq for the continued support to their research work, especially through research projects 308310/2010-1 (A.M.G. de Lima) and 482963/2009-3 (D.A. Rade). The authors express their acknowledgements to the National Institute of Science and Technology of Smart Structures in Engineering, jointly funded by CNPq, CAPES and FAPEMIG.

## References

Austin, E.M., (1999). Variations on modeling of constrained-layer damping treatments. Shock and Vibration Digest 31: 275-280.

Bagley, R.L., Torvik, P.J., (1979). A generalized derivative model for an elastomer damper. Shock and Vibration Bulletin 49: 135-143.

Balmès, E., Germès, S., (2002). Tools for viscoelastic damping treatment design: Application to an automotive floor panel. Proceedings of the 28th International Seminar on Modal Analysis, Leuven, September.

Bathe, K.J., (1996). Finite element procedures. 1. Ed. Prentice Hall.

Bouazzouni, A., Lallement, G., Cogan, S., (1997). Selecting a Ritz basis for reanalysis of frequency response functions of modified structures. Journal of Sound and Vibration 199: 309-322.

de Lima, A.M.G., Stoppa, M.H., Rade, D.A., (2003). Finite element modelling and experimental characterization of beams and plates treated with constraining damping layers. Proceedings of the 17th International Conference of Mechanical Engineering, São Paulo, Brazil.

de Lima, A.M.G., Stoppa, M.H., Rade, D.A., Steffen Jr., V., (2006). Sensitivity analysis of viscoelastic structures. Shock and Vibration 13: 545-558.

de Lima, A.M.G., Kotinda, G.I., Rade, D.A., Steffen Jr., V., Baars, E., (2007). Vibration attenuation in refrigerators using constrained viscoelastic layers: preliminary finite element modeling and experimental assessment. Proceedings of the 22th International Conference on Compressors and their Systems, London, September.

de Lima, A.M.G., Rade, D.A., Lépore-Neto, F.P., (2009). An efficient modeling methodology of structural systems containing viscoelastic dampers based on frequency response function substructuring. Journal of Mechanical System and Signal Processing 23: 1272-1281.

de Lima, A.M.G., da Silva, A.R., Rade, D.A., Bouhaddi, N., (2010). Component mode synthesis combining robust enriched Ritz approach for viscoelastically damped structures. Engineering Structures 32: 1479-1488.

Galucio, A.C., Deü, J.F., Ohayon, R., (2004). Finite element formulation of viscoelastic sandwich beams using fractional derivative operators. Computational Mechanics 33: 282-291.

Golla, D.F., Hughes, P.C., (1985). Dynamics of viscoelastic structures - A time domain finite element formulation. Journal of Applied Mechanics 52: 897-906.

http://www.3m.com/ (Accessed 01 August 2013)

Khatua, T.P., Cheung, Y.K., (1973). Bending and Vibration of Multilayer Sandwich Beams and Plates. International Journal for Numerical Methods in Engineering 6: 11-24.

Kirsch, U., Bogomolni, M., (2007). Nonlinear and dynamic structural analysis using combined approximations. Computers and Structures 85: 566-578.

Lesieutre, G.A., Bianchini, E., (1995). Time domain modeling of linear viscoelasticity using anelastic displacement fields. Journal of Vibration and Acoustics 117: 424-430.
Maia, N.M.M., Montalvão e Silva, J.M., (1997). Theoretical and experimental modal analysis. Research Studies Press LTD., England, p. 468.

McTavish, D., Hughes, P.C., (1993). Modeling of linear viscoelastic space structures. Journal of Vibration and Acoustics 115: 103-115.

Nashif, A.D., Jones, D.I.G., Henderson, J.P., (1985). Vibration damping. John Wiley & Sons, New York.

Rao, M.D., (2001). Recent applications of viscoelastic damping for noise control in automobiles and commercial airplanes. Proceedings of the Emerging Trends in Vibration and Noise Engineering Symposium, Columbus, USA.

Samali, B., Kwok, K.C.S., (1995). Use of viscoelastic dampers in reducing wind and earthquake induced motion of building structures. Engineering Structures 17: 639-654.

Zghal, S., Bouazizi, M.L., Bouhaddi, N., (2014). Reduced-order model for non-linear dynamic analysis of viscoelastic sandwich structures in time domain. MATEC Web of conferences, DOI: 10.1051/matecconf/20141608003.