Finite Element Reduction Strategy for Composite Sandwich Plates with Viscoelastic Layers

Composite materials have been regarded as a convenient strategy in various types of engineering systems such as aeronautical and space structures, as well as architecture and light industry products due to their advantages over the traditional engineering materials, such as their high strength/stiffness relation characteristics and their anti-corrosion properties. This paper is devoted to the finite element modeling of composite laminated structures incorporating viscoelastic materials to the problem of vibration attenuation. However, the typically high dimension of large finite element models of composite structures incorporating viscoelastic materials makes the numerical processes sometimes unfeasible. Within this context, emphasis is placed on a general condensation strategy specially adapted for the case of viscoelastically damped structures, in which a constant (frequency- and temperature-independent) reduction basis to be enriched by static residues associated to the applied loads and the viscoelastic forces is used. After presenting the theoretical foundations, the numerical applications of composite plates treated by viscoelastic materials are addressed, and the main features of the methodology are discussed.


Introduction
In the context of new developments of many types of engineering products mechanical, it is observed that the composite materials have been intensively investigated lately, since they present great advantages over the traditional engineering materials, such as high strength/stiffness relation characteristics and anti-corrosion properties 1 . As a result, the composite materials have been applied in a large number of engineering systems such as automobiles, airplanes, communications satellites, civil engineering structures, etc. However, the most distinguishable feature of the composite structures is that their mechanical behavior depends strongly on a number of factors such as the materials properties and structural configurations, making the numerical modeling a complex task 2 . This is a reason for which in the last decades, a great deal of effort has been devoted to the development of finite element models for characterizing the mechanical behavior of composite materials, accounting for its typical variations of constructions and various orientations possibilities 1 . Much of the knowledge available to date is compiled in the works by Reddy 1 and Berthelot 2 and in some papers such as those by de Lima et al. 3 , Meunier and Shenoi 4 , Chee et al. 5 and Lo et al. 6 .
In engineering applications of composite materials in which dynamic loads are involved, the interest in achieving vibration attenuation becomes of capital importance as vibration amplitudes are directly related to fatigue and, as a result, to structural integrity 3 . Among the various techniques for vibration control which have been devised, the use of viscoelastic materials to be combined with composite structures is an interesting strategy to be investigated. Moreover, for finite element (FE) models of composite viscoelastically damped structures of industrial interest composed by a large number of degrees of freedom (DOFs) the time to compute the exact evaluations during the iterative processes, performed on the full FE matrices, can become prohibitive 3 . Those difficulties motivate the use of a condensation strategy specially adapted for composite viscoelastically damped structures with the intent to reduce the size of the systems 7 .
It is known that fiber-reinforced composite materials present inherent damping mechanisms associated to the viscoelastic behavior of the polymeric matrices and other internal dissipation mechanisms 1,2 . Nonetheless, in a number of cases, such damping may be found to be insufficient to provide the necessary vibration mitigation and must be increased by introducing additional dissipation, which can be done, for example, by applying viscoelastic materials in the form of surface treatments or internal layers 3 . This strategy is considered in this paper.
Among the various theories which have been developed for modeling layered composite structures, the so-named Higher-order Shear Deformation Theory -HSDT proposed by Lo et al. 6 has been chosen in the present study. Despite a more involved analytical framework and an increased number of finite element DOFs, when compared to other simpler theories, the main advantages of the HSDT, which justify such choice, are: (i) it is well adapted to model both thin and thick laminated composite plates; (ii) it accounts for a more complete strain state. As it takes into account transverse shear effects and predicts a parabolic distribution of transverse shear strains, it does not require any correction factor for the distribution of transverse shear strains, such as those required by the First-order Shear Deformation Theory -FSDT. Moreover, it accounts for transverse normal strains.
In the remainder, various theoretical aspects are first presented including: (i) a review of the FE modeling procedure of multilayered plate structures containing unidirectional fiber reinforced and isotropic viscoelastic layers; (ii) the inclusion of the viscoelastic effect into the structural matrices using the complex modulus approach, accounting the frequency-and temperature-dependence of the properties of the viscoelastic material; (iii) the formulation of the constant (frequency-and temperatureindependent) reduction basis specially adapted for viscoelastically damped structures enriched by static residues associated to the external loads and the viscoelastic damping forces. For illustration, the description of two numerical applications of composite plates incorporating viscoelastic materials demonstrates the effectiveness of the viscoelastic damping and the condensation strategy of composite viscoelastically damped structures.

Review of the FE Modeling of Laminated Composite Plates
In this section the formulation of a composite plate finite element is summarized, based on the original developments made by Chee et al. 5 , and implemented by de Lima et al. 3 . Figure 1 depicts a rectangular element formed by k-th unidirectional layer, whose dimensions in directions x and y are denoted by a and b, where z k , h k and θ k indicate, respectively, the thickness coordinate, the thickness and the fiber orientation angle of each layer.
The discretization of the displacement fields within the element is made by using the appropriate interpolation functions through the general relation u(x,y,t) = N(x,y)u (e) (t), where N(x,y) is the matrix formed by the shape interpolation functions, matrices B(x,y,z) are formed by differential operators appearing in the strain-displacement relations 3 , with i = 1 to 8 represents the vector containing the 11 mechanical nodal variables as a function of time. The strain-displacement relations, ε(x,y,z,t) = B(x,y,z)u (e) (t), are used and the resulting strains for elastic layer are generated, where ε xx = ∂u/∂x, ε yy = ∂v/∂y, ε yy = ∂w/∂z, γ xy = ∂u/∂y + ∂v/∂x, γ yz = ∂v/∂z + ∂w/∂y, and γ zx = ∂u/∂z + ∂w/∂x.
Using the displacement and strain interpolations, the expressions of the kinetic and strain energies of the composite plate can be formulated, respectively, as follows 3 : where the elementary mass and stiffness matrices are expressed, respectively, as: x y z dx dy dz where the relation between the displacements at an arbitrary point of the element and vector u(x,y,t) is indicated by matrix A(z). In Equation 2 matrix C (k) (θ k ) represents the orthotropic elastic matrix of the k-th layer, which are constructed according to the Classical Laminate Theory (CLT) as follows: Figure 1. Illustration of the laminated composite plate.
where ( ) k C is the elastic property matrices of the k-th layer, referred to its principal orthotropic axis, and T(θ k ) is the associated rotation matrix 1 .
From the elementary matrices computed for each element of the FE mesh, the global equations of motion are constructed, accounting for the node connectivity, using standard finite element assembling procedures 8 . After assembling, the global equations of motion in the time domain can be written as follows: where ( ) The time domain equations of motion (4) can be used to perform various types of dynamic analysis such as the computation of time domain responses, eigenvalue and frequency response analyses. This is considered later in this paper.

Composite Sandwich Plates with Viscoelastic Layers
The theory presented in the previous section can be easily adapted to the case of sandwich plates containing both unidirectional composite layers and isotropic viscoelastic layers, which is the case considered herein. For these later, constitutive laws must conveniently account for the viscoelastic behavior. Clearly, the unidirectional fiberreinforced layers can also exhibit viscoelastic behavior associated to the inherent behavior of polymeric matrices. Such possibility will be considered in one of the examples considered herein.
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 9 . Various mathematical models have been developed to represent this behavior and have been shown to be particularly suitable to be used in combination with finite element discretization 9 . 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 3,9 . According to the Frequency-Temperature Superposition Principle, also known as Williams, Landell and Ferry (WLF) Principle 10 , 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 by the following expressions: where T indicates an arbitrary value of the temperature, T 0 is a reference value of temperature, ω r = α T (T)ω is the reduced frequency, ω is the actual excitation frequency, and α T (T) is the shift function. Functions G(ω r ) and α T (T) can be obtained from experimental tests for specific viscoelastic materials. In this context, Drake and Soovere 11 suggest analytical expressions for the complex modulus and shift factor the 3M™ ISD112 that is a rubber-like polymer which is provided by the manufacturer in the form of adhesive. Equations 6 represent the complex modulus and the shift factor functions defined in the following temperature and frequency intervals 210 ≤ T ≤ 360 K, 1.0 ≤ ω ≤ 1.0 × 10 6 Hz, respectively, where T 0 = 290 K, for the 3M™ ISD112 viscoelastic material 12 , as given by those authors. The parameters appearing in the following expressions are presented in Table 1.
From the reduced temperature nomogram generated by the computation of expressions (6), the designer can obtain the complex modulus and the loss factor at any given temperature into a frequency band of interest, as illustrated in Figure 2 for the particular material considered above 13 .
According to the Elastic-Viscoelastic Correspondence Principle 10 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 by the underlying theory. This step has been accomplished in the preceding section. Then, the material moduli are modified to account for the viscoelastic behavior (according to the model expressed by Equation 6). Clearly, this approach leads to frequency-and temperature-dependent FE stiffnesses matrices, which are expressed, after a proper adaptation of Equation 2b, as follows: x y z T x y z dx dy dz At this step, one of the moduli can be factored-out of the viscoelastic stiffness matrix that represent its contribution to the stiffness matrix of the sandwich plate as follows: is frequency-and temperature-independent matrix, which is combined with the stiffness matrix represented by the Equation 2b, to produce the following complex global stiffness matrix of the composite sandwich plate with viscoelastic material: where K e and K v (ω,T) are, respectively, the contributions of the purely elastic and viscoelastic layers to the global stiffness matrix. Neglecting other forms of damping, after Fourier transforming of Equation 4, the global finite element equations of motion in the frequency domain of the composite sandwich plate can be expressed as follows: where N is the total number of DOFs. M ∈ R N × N is the mass matrix and K e ∈ R N × N and ( ) are the stiffness matrices corresponding to the elastic and viscoelastic layers, respectively. Q(ω,T) ∈ R N and F(ω) ∈ R N are, respectively, the vectors of the amplitudes of the harmonic generalized displacements and external loads. The so-named receptance or frequency response function matrix (FRFs), the components of which are complex functions of the excitation frequency that establish linear relations between the amplitudes of the harmonic responses and the amplitudes of the excitation forces and moments is expressed as: , -

Reduction of Composite Structures with Viscoelastic Material
In the case of complex composite structures of industrial interest, FE models are usually constituted by a large number of DOFs. In such cases, it becomes practically impossible to compute the FRFs directly from Equation 11, owing to the prohibitive computation times and storage memory required. This fact motivates the use of model condensation 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 expression (10), can be approached by projections on a reduced vector basis as:

Q(ω,T) = TQ(ω,T)
where T ∈ C N × NR is the transformation matrix formed column-wise by a vector basis, Q(ω,T) ∈C NR are generalized coordinates, and NR << 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 reduced model provides a reasonable predictive capacity into a frequency bandwidth. It must be emphasize that the 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 Equations 10 and 12, the FRF matrix (11) can be rewritten as follows: (13) For models containing viscous or structural damping, it is relatively common to use a constant projection basis formed by the eigenvectors of the associated conservative structure, as the mass and stiffness matrices are invariant 7,9 . 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. 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 by Balmès and Germès 14 , 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 η 0 = 0. On the other hand, for high frequencies, the extrapolation gives the complex values G∞ and η∞. The tangent stiffness matrix can thus be calculated as follows 13 :

H T K T T K T T MT
The nominal basis can be obtained by the resolution of the eigenvalue problem: This basis is enriched by introducing the residues formed by the static displacements associated to external forces, that must be further completed by the residual vectors associated to the viscoelastic damping forces: where matrix b ∈ R N is a Boolean matrix which enable to select, among the DOFs, those in which the excitation forces are applied. These residuals are interpreted as the columns of the flexibility matrix of the associated undamped system, related to the coordinates of application of two types of forces to it: the external excitation forces and the viscoelastic damping forces. These latter can be better understood by examining Equation 10, 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: 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 7 . Also, experience has demonstrated that the nominal basis (18) 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 7,14 .

Condensation of composite structures with inherent structural damping
In this section, numerical tests were performed using the FE model of a square composite plate without viscoelastic material depicted in Figure 3   It is assumed that the composite material present inherent hysteretic damping, represented by complex, frequencyindependent moduli of the form ( ) , in which a single value of the loss factor η mn = 0.001 is adopted for all the moduli 3 .
The interest is to evaluate the nominal enriched basis of reduction by using the static residues associated to the applied external loads. The computations consisted in obtaining the driving point FRFs  Figure 4 show the amplitudes of the FRFs computed by using the nominal basis, as compared to the amplitudes of the dynamic response computed by using a reference basis formed by a far larger number of eigenvectors (400) and one residual vector associated to the applied external load. Also, the FRF amplitudes have been computed by using a convenient reference factor through the relation It can be clearly seen that the accuracy is continuously improved upon enrichment of the reduction basis by the inclusion of the static residues associated to the external loading to form the nominal basis T 02 . Thus, the use of first order residues associated with external forces are sufficient to represent with accuracy the dynamic behavior of composite laminated structures with inherent structural damping.

Condensation of composite sandwich structures with viscoelastic layers
The interest now is to evaluate the accuracy of the nominal basis further enriched to account for viscoelastic damping forces for composite structures incorporating viscoelastic materials. It is considered a sandwich rectangular plate composed by 4 unidirectional fiberreinforced composite layers and a polymeric core made of viscoelastic material 3M™ ISD112, with modulus given by Equations 6 and mass density ρ = 950 kg.m -3 . The FE discretization mesh, the geometrical characteristics and the boundary conditions of the composite sandwich structure considered here is the same as those of the plate illustrated in Section 5.1 (see Figure 3), with the exception that the middle layer consists of the viscoelastic material with thickness h v = L x /128 m. The upper and the bottom fiber-reinforced layers have the same thickness and material properties as those considered in the previous section.
To verify the condensation, one considers the following nominal basis: ) after SVD filtering). The residues R v 0 were computed based on the largest singular values, for which the relation σ 1 /σ i ≤ 1 × 10 5 , for i = 1 to 10 holds. Figure 5 shows the FRF amplitudes computed by using the three nominal bases as compared to the amplitudes of the dynamic response of the reference system evaluated for a nominal temperature value of the viscoelastic material of 25 °C. It is shown that the enrichment of the classical basis of reduction associated to the external loads and the viscoelastic damping effects improves the results in the frequency band of analysis, leading to conclude that approximations based only on the use of static residues due to the external loads might not be accurate enough to predict the dynamic behavior of the viscoelastically damped composite structures.

Concluding Remarks
A condensation procedure intended to be used for dealing with viscoelastically damped composite laminated structures was suggested and evaluated. The original aspect of the procedure lies in the adaptation of the concept of condensation, initially developed for undamped structures, for composite structures containing viscoelastic materials to the problem of vibration attenuation. This approach allows design procedures of real-word complex engineering structures made of composite materials with viscoelastic damping, for example in the context of an optimization and/or model updating processes, without the necessity to compute the dynamic responses of the complete system. Thus, this condensation strategy has been proposed aiming at alleviating the computational costs involved in the analyses based on large-scale FE models of composite structures containing a large number of DOFs.
The academic examples were used to illustrate the efficiency of the condensation procedure, mainly in terms of the drastic reduction of the number of DOFs, which demonstrates that the suggested technique is adapted to more complex viscoelastically damped composite systems, and being a very useful tool for the design, analysis, structural optimization and/or model updating processes.