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.


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-1 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.

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 1 u and 1 v , respectively), the in-plane displacements of the middle plane of the constraining layer in directions x and y (denoted by 3 u and 3 v , respectively), the transverse displacements, w , and the cross-section rotations about x and y , denoted by x  and 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, x  and y  , are the same for the elastic layers; (vi) the 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, x y N is the matrix formed by the shape functions, and , are used and the resulting strains for the elastic layers, k  ε , and for the viscoelastic core, 2  ε 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: is the isotropic elastic material properties matrix, and 2 v    C C contains the frequency-and temperature-dependent material properties for the viscoelastic core. Matrix

 
x, y, z B 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 dimen-sions of the plate element in directions x and y, and k h and 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, , 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: 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: (1) where the stiffness matrices,   give, respectively, the contributions of the purely elastic and viscoelastic parts of the damped structure, and 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.

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: where 0 E and E  are, respectively, the relaxed or static modulus, and non-relaxed or highfrequency limit value of the modulus,  is the relaxation time, and  represents the fractional order of the time derivative   0 1    . One important aspect regarding the use of the following complex modulus function, , obtained by applying the Fourier transform to Eq. (4) is the identification of the four parameters ( 0 E , E  ,  ,  ) from experimental data-sheets provided by the manufactures 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. 6 Also, it must be emphasize that Eq. (4) was obtained by introducing the following internal variable as a strain function,       , in the standard one-dimensional stress state constitutive equation for linear viscoelastic materials. As a result, Eq. (4) contains only one fractional derivative term, , instead of two presented in the classical constitutive equation, and can be approximated by the following Grünwald definition, , noting that 1 1 A  , to generate the following discretized form of the anelastic strain: is the time step, p n n  is the number of discretization points, and 1 j A  represents the Grünwald coefficients given by the recurrence formula, appearing in equations of motion (3) can be expressed as, where for a general stress-strain relation, the vector   , , , x y z t σ can be written in terms of the anelastic strain (5) taking into account the strain relation, , as follows: By combining Eq. (6) with the elastic and anelastic strain-displacement relations, ε u x, y, z,t x, y, z , the resulting modified elementary equation of motion takes the following form: where 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: is a load vector, which depends on the anelastic displacements,   t u , 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][6][7][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.

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: is the transformation matrix formed column-wise by vectors,   NR t C  u is the vector of generalized coordinates, and NR 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: 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, 0 E , as follows: Hence, the nominal basis can be obtained by the resolution of the eigenvalue problem: The basis, 0 ϕ , is further enriched by introducing the following residues formed by the static displacements associated, respectively, to viscoelastic damping forces and external excitations: 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: (13)

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: 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: , of the linear behavior, and can be view 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: where 0 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, i  ϕ is the perturbation due to the local nonlinearities which can be interpreted as an external load applied on the linear model, 0 i  ϕ . Thus, a set of static responses can be generated: 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. count the reduced matrices, M and 0  K ; (iv) in the next step, a new set of reduced velocities,

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 baseplate 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×10 9 N/m 2 ; mass density   2700 kg/m 3 ; Poisson ratio,   0.3. For the 3M™ ISD112,   950kg/m 3 ,   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,   , where 0 F  10 N is the amplitude of the excitation and 0 f is the excitation frequency, and to a transverse triangular impulse loading applied to point I indicated in Fig. 4  . For the harmonic loading, the excitation frequency ( 0 f  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, 0 E . Also, it was considered a proportional damping computed by the following relation, 0.005   D 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: 1 (19 eigenvectors, plus one residual vector computed by the Eq. (13.b)); (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 0 v R were computed based on the largest singular values, for which the relation, max i    1x10 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 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. 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.

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.

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 nonlinear 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 that the efficiency of the proposed timecondensation 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.