Dynamic Analysis of Three-Layer Sandwich Beams with Thick Viscoelastic Damping Core for Finite Element Applications

This paper presents an analysis of the dynamic behaviour of constrained layer damping (CLD) beams with thick viscoelastic layer. Ahomogenisedmodel for the flexural stiffness is formulated usingReddy-Bickford’s quadratic shear in each layer, and it is compared with Ross-Kerwin-Ungar (RKU) classical model, which considers a uniform shear deformation for the viscoelastic core. In order to analyse the efficiency of both models, a numerical application is accomplished and the provided results are compared with those of a 2Dmodel using finite elements, which considers extensional and shear stress and longitudinal, transverse, and rotational inertias. The intermediate viscoelastic material is characterised by a fractional derivative model, with a frequency dependent complex modulus. Eigenvalues and eigenvectors are obtained from an iterative method avoiding the computational problems derived from the frequency dependence of the stiffness matrices. Also, frequency response functions are calculated. The results show that the new model provides better accuracy than the RKU one as the thickness of the core layer increases. In conclusion, a new model has been developed, being able to reproduce the mechanical behaviour of thick CLD beams, reducing storage needs and computational time compared with a 2D model, and improving the results from the RKU model.


Introduction
In the last years many studies have been presented concerning the structural vibration reduction making use of passive damping control techniques by means of surface treatments with viscoelastic materials.A survey of different subjects on viscoelastic treatments can be found in [1].This kind of vibration control technique is largely used nowadays for several industrial applications, such as aeronautical and automotive components.The free layer damping (FLD) and constrained layer damping (CLD) technologies are two of these viscoelastic surface treatments, consisting of adding a damping viscoelastic layer to the structural system.Specifically, FLD consists of adding that viscoelastic layer on a vibrating metallic base, and the configuration can be analyzed as the flexural behavior of a two-layer beam.In this context, Oberst and Frankenfeld's model [2] is traditionally used.This model assumes that the layers are thin and the shear effect is negligible.In one of the authors' works [3] a model was presented in which shear effects are taken into account, improving the results given by Oberst and Frankenfeld's model for thick layers.On the other hand, CLD technologies are based on a viscoelastic core working under shear deformation.In this configuration, this damping viscoelastic material is in the core of a three-layer sandwich structure.The other two layers are vibrating structural elements, usually made of metallic materials subjected to flexural moment and axial and shear loads.The effect of this shear load is more important as the thickness of the layers increases.Compared with other damping techniques, this viscoelastic treatment presents the disadvantage of adding mass into the original system but offers a better damping-weight ratio than the free layer damping configuration [4,5].
In Figure 1 the cross section of a three-layer composed beam is represented, the thickness of the base, viscoelastic, and constraining layers being  1 ,  2 , and  3 , respectively, and  being the width of the beam.Although in CLD beams  2 ≪  1 , in this work three-layer sandwich beams with any value for  2 are considered.The properties for these elastic materials are Young's moduli  1 and  3 , Poisson's coefficients ] 1 and ] 3 , and densities  1 and  3 , whereas those of the viscoelastic core are  * 2 , ] 2 , and  2 , respectively.If the materials are considered isotropic, the shear moduli  1 ,  * 2 , and  3 are related to the extensional ones by , and  3 =  3 /2(1 + ] 3 ) for the three layers.
The complex nature of the modulus of viscoelastic materials, represented by the asterisk as (⋅) * , is due to their ability to dissipate mechanical energy.A complex relationship between stress and strain allows representing the hysteretic behaviour of this kind of materials [6].The dissipative property may be described by the loss factor  2 : where   2 and   2 are the real and imaginary components of the complex modulus  * 2 =   2 +   2 , which are known as the storage and loss modulus, respectively, the former being also frequently represented by  2 .For polymers, both storage and loss moduli and consequently loss factor vary with temperature and frequency in terms of three different behaviours: rubbery, vitreous, and transition (see, e.g., [7]).In this sense, Jones [8] reviews experimental data for complex modulus of some typical viscoelastic materials, including elastomers, adhesives, and specially compounded materials.The experimental characterisation of low stiffness damping materials, which are not appropriate to prepare selfsupporting test specimens, may be achieved by means of the ASTM E 756-05 standard [9], wherefrom the nomenclature has been taken.
The classical RKU [10] model is one of the first models developed in order to calculate the flexural stiffness for sandwich beams with viscoelastic core and is one of the most used models nowadays, although it may lose accuracy in several cases, for example, when rotational or extensional inertias are not negligible and specifically when a thick viscoelastic core in a CLD beam is studied.This is due to the assumptions which this model takes.
(i) The shear deformation in the elastic base and constraining layers is considered to be negligible.
(ii) The longitudinal direct stress in the viscoelastic layer is negligible.(iii) The shear stress in the viscoelastic layer is assumed to be uniform.(iv) The plane cross section in each layer remains plane after deformation.
There are some other models beyond these restrictions such as [11][12][13][14][15][16][17][18] or more recently [19,20] but they are not so computationally efficient for finite elements applications.In fact, RKU is used in most engineering applications basically due to its very easy computational implementation.Also, RKU strictly applies to simply supported beams, although it is much more generally used, for other boundary conditions, such as the fixed support.We will refer to the results obtained as RKU results, although RKU theory for simply supported beams has been used in this work for a cantilever beam, in the context of a finite element analysis with displacement fields restricted to those of the RKU theory.
In short, following a similar approach as for the study of the FLD configuration for thick beams previously mentioned [3], this work is aimed at developing a new model improving the accuracy of the RKU model for three-layer sandwich beams with a thick viscoelastic core but maintaining its computational benefit.Thus, an equivalent complex flexural stiffness is derived considering quadratic shear stress based on Reddy-Bickford's theory.In order to prove the improvement achieved by the new model, a numerical application for a cantilever beam is presented comparing the solutions provided by three different finite element models: a 2D model (whose results are considered to be the reference ones) and two 1D models, based on the RKU theory and the new one, respectively.The damping material is characterised by means of a fractional derivative model involving the frequency dependence of the complex modulus, which implies important disadvantages for the dynamic analysis.Thus, the extraction of the eigenvalues and eigenvectors is carried out by a simple and effective iterative algorithm [21].Finally, the dynamic response of the three models is compared in terms of the frequency response function.
As a result, the new model is able to reproduce the mechanical behaviour of three-layer sandwich beams, reducing storage needs and computational time compared with a 2D model and improving the results provided by the classical RKU model.

Homogenised Model for a CLD Beam
Next the theoretical study of a three-layer sandwich beam is presented, in which quadratic shear stress is taken into account.An equivalent flexural stiffness will be deduced for pinned-pinned beams.In order to simplify the notation, it is assumed that the behaviour of all materials is linear and elastic (i.e., all the magnitudes are real) and the complex character of the modulus in the viscoelastic core will be taken into account in the numerical examples of the subsequent sections.This substitution of a complex modulus into an elastic solution is usually known as the "correspondence principle" (see any standard text on viscoelasticity, e.g., [22]).
The equivalent flexural stiffness  eq may be obtained by the addition of the individual contribution of each layer: where  1 =  1  1 ,  2 =  2  2 , and  3 =  3  3 are the complex flexural stiffness of the three layers and  1 ,  2 , and  3 are the complex cross-sectional second order moments, computed with respect to the neutral axis, given by respectively, in which the complex position of the neutral axis is represented by Following the same methodology as in [3] and decoupling the transverse displacement of any cross section in a term due to the flexural moment and in another term derived by the shearing force (see any book of strength of materials for details, e.g., [23]), the equivalent shear stiffness of the cross section can be found to be decomposed as For the geometry represented in Figure 1, the stiffness of  1 ,  2 , and  3 of the individual layers satisfies respectively, where respectively.In these equations,  eq is the flexural stiffness given by (2).By solving these integrals, it yields respectively, where For base and constraining layers made of metallic materials, which are much thinner and more rigid than the viscoelastic layer, the terms 1/ 1 and 1/ 3 may be neglected with respect to the coefficient of the polymeric layer 1/ 2 .
It can be pointed out that if the three layers were composed of the same material, ( 14) provides the well-known result for homogeneous rectangular sections: where  is the total cross-sectional area.By considering the shear coefficient  eq , the flexural field equation in free vibration is given by Timoshenko's formula [24]: where  ℓ is the mass per unit length.In (18) the rotational inertia has been neglected.When the shear stiffness tends to infinity, (18) degenerates on the well-known Euler-Bernoulli field equation.Following the same methodology as in [3], the equivalent flexural stiffness   , considering shear effects and satisfying the homogenised Euler-Bernoulli field equation, can be found to be The function () takes into account the shear effects and is given by In order to compare RKU and the new model, it can be noted that, for the static values,  → 0, both models give the same result, the equivalent  eq .However, as a difference, when frequency tends to infinity  → ∞, the stiffness of the present model tends to zero whereas that of the RKU one tends to a finite value.These two models have another common property when shear stiffness  eq tends to infinity, that is, when shear deformations are negligible.In this case, the function () tends to zero and the equivalent flexural stiffness   tends to the classic  eq .

Dynamic Analysis of a Three-Layer Sandwich Beam Using Finite Element Procedures
3.1.Problem Definition.In this section the harmonic analysis of a three-layer sandwich beam in a CLD configuration will be completed using finite element procedure techniques.Three different thicknesses of the viscoelastic core layer will be studied, so as to evaluate the accuracy of the homogenized models compared to a 2D model, whose solution will be considered as exact, considering the nonexistence of experimental results.The length of the beam is ℓ = 200 mm, the width is  = 20 mm, the thickness of the base and constraining metallic layers is  1 =  3 = 1 mm, and for the viscoelastic layer  2 = 1, 5, and 10 mm is chosen.
The properties of the materials are taken from the experimental characterisation effectuated by Cortés and Elejabarrieta [25] on AISI T 316 L stainless steel laminated sheet and on Soundown Vibration Damping Tile material [26].Indeed, Young's modulus and density of the elastic material are  1 = 176.2× 10 9 Pa and  1 = 7782 kg/m 3 , respectively, and density of the damping material is  2 = 1423 kg/m 3 .Poisson's coefficient ] 1 = ] 3 = 0.3 is chosen for the elastic materials, and ] 2 = 0.45 is chosen for the viscoelastic material.The experimental data of the storage modulus and loss factor for the viscoelastic core layer were fitted to a fourparameter fractional model [27,28], given by where   and   represent the relaxed and unrelaxed modulus, respectively;  is the relaxation time and  is the fractional parameter.The parameter values are summarised in Table 1.The dynamic behaviour of the three-layer sandwich beam in a CLD configuration is studied on the basis of three different finite element models.The first of them is a 2D model discretised in bilinear quadrilateral elements with four nodes under plane-stress assumption (see, e.g., [29][30][31] for details about finite element formulations).All the three layers are modelled with 4 elements along thickness to assure the continuous evolution of the shear stress and with 60 elements along the length (see Figure 2) in order to obtain the first three eigenvalues accurately enough: it has been checked that this number of finite elements is enough to get convergence for any of the results shown in the tables.
The consistent mass matrix and the stiffness matrix obtained using reduced integration with Kosloff and Frazier [32] hourglass control are summarised as follows.
Mass M and Stiffness K Matrices for the 2D Model Quadrilateral Finite Elements.See Figure 3.
Stiffness Matrix K Obtained by Reduced Integration with Hourglass Control.Consider Extensional and shear strain and transverse, extensional, and rotational inertias are considered in this 2D finite element model which is then able to reproduce the dynamical behaviour of the three-layer sandwich beam with any thickness of the viscoelastic core.
The two other models are 1D beam models whose th mass M  and stiffness K  matrices are respectively, where   is the length of the th finite element.The complex flexural stiffness  of ( 25) is given by ( 2)  =  * eq for the RKU model and by ( 19)  =  *  for the new thick beam model.The discretisation is also made with 60 finite elements along span.

Extraction of Eigenvalues.
The equation from which the complex eigenvalues of the system under study can be obtained is where  *  and  *  are the complex eigenvalue and eigenvector of the th mode, respectively, M is the mass matrix, and K * is the complex stiffness matrix, which is dependent on frequency.This   is the real part of the square root of the complex eigenvalue  *  : which induces a nonlinearity into the eigenproblem.There are several methods such as those of Lanczos [33] or Arnoldi [34], which use iterative procedures involving important computational time.In order to decrease this computational effort, Cortés and Elejabarrieta [21] developed an iterative procedure that approximates in a simple and accurate way the complex eigenpair.This method begins by considering the static stiffness matrix K * (0) in (26) yielding and then the undamped eigensolutions  ,0 and  ,0 can be obtained.Then, Nelson's method [35] is used in order to calculate the eigenvector derivative    .From this eigenvector derivative and taking into account the variation of the complex stiffness for the obtained eigenfrequency, and by means of Taylor's series approach, a complex finite increment Δ *  of the eigenvector is obtained.The complex eigenvector can be approximated with with which the complex eigenvalue  *  is estimated according to where (⋅)  denotes the Hermitian transpose operator, that is, the complex conjugate transposition.Equations ( 29)-( 31) can be iterated making use of the new eigenfrequency   given by (27), in order to obtain the desired convergence tolerance.As a main difference with other iterative methods, this one presents the advantage of solving only once the undamped eigenproblem, and the iterations are carried out on the derivatives, reducing computational resources.If damping in the system is very large, the accuracy of the method can be improved by means of the incremental approach of the method, as seen in [21].
Making use of this new method, with the corresponding incremental approach, the first three modal natural frequencies   derived from (27) and loss factor derived from can be computed.The corresponding results for the three thicknesses and for the three models under study are shown in Tables 2-4.
It can be pointed out that the results for the natural frequency   are practically the same for the thinnest beam (see Table 2), and the differences between the present model and the RKU one are more important as the thickness of the viscoelastic layer increases, which is an expected behaviour; the reason is that the shear contribution is more important for larger thickness, and the present model considers a more realistic shear stress distribution.Also, the most important differences take place at higher order modes.This is because, at higher frequencies, the shear effects acquire more importance, and, as previously mentioned, the present model takes into account shear effects in a more effective way.
Specifically, for the third mode of the beam with a viscoelastic layer thickness equal to 5 mm (see Table 3), the present model improves the RKU result in a 6.8% (from 10.5% down to 3.7%).This improvement is even more important for the beam with 10 mm of viscoelastic layer (see Table 4), where the difference between the present model and the RKU one with respect to the 2D model goes up to 16.5% (from 21.5% to 5.0%).
As for the results of the modal loss factor   , an erratic behaviour in both RKU and thick beam models can be noted.It should be highlighted that this parameter cannot be directly compared, because the damping of the viscoelastic material represented by the loss factor   depends on frequency according to (21), and the natural frequencies for the models are not the same.Instead of modal loss factor   , the amplitudes of the resonance peaks will be compared in the next section.

Analysis of the Frequency Response
Function.The matrix system for a cantilever beam is given by where F * and U * are the nodal vectors of the amplitude of the force and the displacement, respectively.The frequency dependence of the complex stiffness matrix K * involves the fact that the modal superposition cannot be applied, and therefore (33) must be solved for U *  at each desired frequency   .Figure 4 represents the frequency response up to 3 kHz of the free-end of the three-layer sandwich beam in all the three models, when a harmonic unitary force  = 1 × exp()  is applied on the right-hand side of our system.
Figure 4(a) indicates that the responses provided by the three models are practically the same, as expected for such a thin beam in which the shear is not so important.On the contrary, Figures 4(b) and 4(c) illustrate a displacement of the curves to the right with respect to that given by the 2D model.This is due to the overestimation of their natural frequencies, as it was mentioned in the previous section.The differences are more patent as the frequency is larger and as the thickness is bigger, mostly for the RKU model.
Regarding the amplitudes of the resonance peaks, Tables 5-7 show how the model presented in this paper improves the results of the RKU model, the improvement being better as the thickness of the viscoelastic layer increases.For instance, for the first mode of the thickest beam, the difference between the errors in RKU and present models goes down from 7.3% to 3.4%, which is a significant difference taking into account the logarithmic scale.
Since a quadratic distribution of the shear stress is representative of the state of stress in beams, the differences between the new model and the 2D model are mainly due to the inertias.

Conclusions
In this paper a formulation for three-layer sandwich beams has been presented, which takes into account the deflection due to shearing forces.For this, quadratic distribution of shear stress was considered along thickness.In order to test the performance of the presented model, a dynamical analysis has been carried out on a CLD beam, in which the viscoelastic material has been modelled by means of a fractional derivative model.Consequently, the natural frequencies and resonance peak amplitudes provided by the new model and by the classical RKU model have been compared with those given by a 2D finite element model.As a conclusion, both models present similar results for a thin core layer, but the results are significantly improved as the thickness of the core layer increases.This formulation is very simple to be implemented in finite element analysis, presenting the same complexity level as the RKU one.This represents an important advantage with respect to other formulations.
In short, a new model for finite element calculations in three-layer sandwich beams with viscoelastic core layer has been developed, improving the results of the RKU one, because the former considers shearing force deflection, without enlarging computational efforts.

Figure 1 :
Figure 1: Cross section of a three-layer sandwich beam with a viscoelastic damping layer.

Figure 2 :
Figure 2: Finite element model for the 2D dynamical analysis of a sandwich beam.
Aspect ratio:  = b/a Material properties: E,  and

Table 1 :
Parameters of the fractional derivative model.

Table 2 :
Modal properties of the sandwich cantilever beam with  2 = 1 mm.

Table 3 :
Modal properties of the sandwich cantilever beam with  2 = 5 mm.

Table 4 :
Modal properties of the sandwich cantilever beam with  2 = 10 mm.

Table 5 :
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thickness  2 = 1 mm.

Table 6 :
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thickness  2 = 5 mm.

Table 7 :
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thickness  2 = 10 mm.