General Homogenised Formulation for Thick Viscoelastic Layered Structures for Finite Element Applications

Viscoelastic layered surface treatments are widely used for passive control of vibration and noise, especially in passenger vehicles and buildings. When the viscoelastic layer is thick, the structural models must account for shear effects. In this work, a homogenised formulation for thick N-layered viscoelastic structures for finite element applications is presented, which allows for avoiding computationally expensive models based on solids. This is achieved by substituting the flexural stiffness in the governing thin beam or plate equation by a frequency dependent equivalent flexural stiffness that takes shear and the properties of the different layers into account. The formulation is applied to Free Layer Damping (FLD) and Constrained Layer Damping (CLD) beams and plates and its ability to accurately compute the eigenpairs and dynamic response is tested by implementing it in a finite element model and comparing the obtained results to those given by the standard for the application—Oberst for the FLD case and RKU for the CLD one—and to a solid model, which is used as reference. For the cases studied, the homogenised formulation is nearly as precise as the model based on solids, but requires less computational effort, and provides better results than the standard model.


Introduction
When a structure is subjected to vibration and the emitted sound must be reduced, surface treatments are a handy solution. This kind of passive control is specially used in appliances, building, and human transportation [1][2][3]. In this last case, apart from reducing structure-borne sound [4], this technology serves the purpose of increasing passenger comfort.
Surface treatments are layered slender structures that can be divided into two groups: free layer damping (FLD) and constrained layer damping (CLD). The first consists of a layer of a base material, usually metallic, to which a layer of viscoelastic material is glued either in the form of damping tiles [5] or coating [4]. This way, when the treated structure deforms in bending, the viscoelastic layer suffers an extensional effort and dissipates energy in the form of heat. CLD surface treatments, on the contrary, have an additional constraining layer, usually metallic as well that forces the viscoelastic core to deform in shear so that energy is dissipated [1]. Even if the CLD configuration is more effective to reduce vibration, FLD treatments continue to be applied due to their easy and relatively economic implementation [6].
In order to analyse the dynamic behaviour of such layered structures, several strategies have been proposed, ranging from analytical to numerical methods. The former approach can be followed when the RKU for those with a CLD configuration. For every case, the dynamic behaviour of the viscoelastic layer is represented by a four parameter fractional derivative model.

General Formulation for Viscoelastic Layered Surface Treatments
In order to develop a general formulation, the N-layered beam shown in Figure 1 is considered. The process goes as follows: first, an equivalent modulus that takes into account the contribution of the different materials to the overall stiffness of the beam is computed; then, the effect of shear is introduced by a frequency dependent term that modifies this equivalent modulus. The formulation for layered beams will be presented and, afterwards, it is extended to plates.  Following the classical theory for layered beams, the total flexural stiffness B eq for a N-layered beam can be expressed as the sum of the flexural stiffness of its layers as where E i stands for the elastic modulus of the material of the ith layer and I i for the cross-sectional inertia moment computed from the neutral axis. This inertia term is obtained for the ith layer from the integral where b is the width of the beam that is supposed to be constant, and defining the integration limits . . .
for every layer, where h n stands for the position of the neutral fibre that is located in If the small deformation hypothesis holds, it is possible to uncouple the transverse displacement due to bending v M (x, t) and to shear The first term, as it is originated by the bending moment M(x, t), satisfies where B eq is the equivalent flexural stiffness from (1). Likewise, the second term is related to the shear force Q(x, t) by where the total shear stiffness K eq is obtained by combining the shear stiffnesses of all layers If a second order shear stress law is considered, the shear stiffness for the ith layer follows where G i is the shear modulus of the material of the ith layer and where S represents the area above a point P(x, y) (see a book on mechanics of materials for the details, e.g., [25]). Taking into account the integration limits defined in (3), Γ i (y) has the value for the bottom layer, for the middle layers 2 to N-1 and for the top layer N.
The equivalent flexural and shear stiffnesses B eq and K eq are then introduced in the Timoshenko thick beam formulation [26] in order to account for the effect of the different materials of the layers, resulting in if the terms related to the rotational inertia of the cross section are ignored, and where q(x, t) stands for the distributed force per unit length applied to the beam and ρ L = ∑ N i=1 ρ i A i being the mass per unit length and ρ i and A i the density and area of the cross section of the ith layer.
If harmonic response both in space and time is assumed, the transverse displacement v(x, t) takes where A, β and ω are the amplitude, wave number and natural frequency. Introducing (15) in (14) and after some manipulation, an equation equivalent to the Euler-Bernouilli thin beam formula [27] but which includes the effect of shear can be obtained whereB(t) is the time domain counterpart of the flexural stiffness modulus that, in the frequency domain, takes the form in which Two considerations should be made regarding this final expression: the first is that, if the modulus of the material of any of the layers is complex, the formulation is still valid, the flexural stiffness terms becoming complex; the second is that, even if the dependence on frequency could be thought as a unnecessary addition of complexity, most of the time, the viscoelastic materials used in these kinds of applications already show frequency dependent properties [28].
It should also be noted that the presented formulation only accounts for the bending behaviour of the structure under study and that is only valid for the dynamic regime i.e., it cannot be used to compute the static response of thick beams in which the shear is relevant. Considering that the domain of application is the reduction of noise and vibration, this does not pose a problem.
In order to extend the formulation to plates, the mass per unit length ρ L should be substituted by the mass per unit area where ρ i and H i stand for the density and thickness of the ith layer. In addition, when dealing with bending, the equivalent flexural stiffness B eq should be substituted by its counterpart for plates E i and ν i being the elastic modulus and Poisson's ratio of the ith layer. The integration limits are the same ones defined for beams in (3). The shear formula still holds taking into account that, for plates, the shear force is considered per unit width and, then, b = 1 in (9) and (11)- (13). This way, the frequency dependent flexural stiffness is given by where ϕ(ω) follows the expression Considering also for plates y the vertical axis and v the transverse displacement in order to avoid the confusion between the typically used w and the natural frequency ω, the Kirchhoff-Love thin plate equation [27] becomes where p(x, z, t) is the distributed force per unit area applied to the plate andD(t) stands for the flexural stiffness in the time domain and introduces the effect of shear into the formulation. The same considerations already mentioned for the beams apply here: the presented procedure is only applicable when working in the frequency domain.

Dynamic Analysis
The presented homogenised formulation is used to compute, first, the natural frequencies and mode shapes of layered beams and plates and, later, their response in the frequency domain. To this aim, the finite element method is used. As the formulation is directly based on the thin beam (or plate) equation, standard elements can be used, the single difference being that the stiffness matrix is, in this case, frequency dependent.
This fact does not add much complexity taking into account that, due to its definition, the formulation allows for computing the frequency dependent stiffness from the static stiffness matrix K(0) as for beams and for plates. This implies that an external program can be used to produce the mesh for the structure, which drastically simplifies the computation process. It is also worth mentioning that, in a general case, the frequency dependent stiffness matrix can be complex-this is the reason why it is marked with an asterisk, in a fashion that will be followed all along the work. This said, the process for computing, on the one hand, natural frequencies and mode shapes and, on the other, the dynamic response is referred.

Natural Frequencies and Mode Shapes
If the finite element method is applied either to (16) or (23) when the structure is not subjected to an external force, the obtainment of the response becomes the eigenvalue and eigenvector problem where M is the mass matrix, ω r = Re λ * r is the natural frequency for mode rth and λ * r and φ * r are the corresponding eigenvalue and eigenvector (see a book dealing with finite element formulation for the details, for example, [29]). As the stiffness matrix is frequency dependent, an iterative method is needed to compute the eigenpairs of the system. For each mode, the computation starts considering the static stiffness matrix K(0) and computing the associated natural frequency. Then, the algorithm shown in Figure 2 is followed, obtaining the natural frequency from the updated stiffness matrix and comparing it to the one computed in the previous iteration. When the difference in natural frequency for two subsequent iterations is less than a given tolerance, the process stops. It should be noted that the iterative method can be expensive in terms of computational resources. In such case, approximate techniques like the one presented in [30] can be applied to reduce the computational cost.
In the case under study, apart from the static case, the eigenvalues of the system are complex due to the fractional model used to represent the behaviour of the viscoelastic material. If the viscoelastic layer covers the whole structure, the eigenvectors are real, due to the proportionality of damping, but, in a general case, both eigenvalues and eigenvectors could be complex.
Apart from the natural frequencies and mode shapes, the loss factor is also used to compare the performance of the different models. For the rth mode, it follows that

Dynamic Response
The dynamic response of the structure is obtained by solving the motion equation in the frequency domain for the response vector v * (ω). For the computations, a uniform impulse pressure (1 Pa) on the upper surface of the structure is considered as the force vector F(ω) and, for the comparison between models, the RMS value of the transverse displacement v * (ω) on that same surface is computed as where R is the number of nodes of the model and v * i (ω) the response in the frequency domain of the ith node in the transverse direction.

Case Studies
The presented methodology is then tested in Matlab (R2019a, The MathWorks Inc., Natick, MA, USA) for both beam and plates, for the FLD and CLD configurations in each case (The code developed can be found as Supplementary Materials). For the sake of brevity, only the results for the simply supported case are reported; cantilever FLD and CLD beams are studied in [21,22] and the behaviour of FLD and CLD plates with different boundary conditions is analysed in [23,24].
The computations are done considering the density of the steel layers ρ 1 = 7782 kg/m 3 and their Young's modulus E 1 = 176.24 GPa. The viscoelastic layer has a density of ρ 2 = 1423 kg/m 3 and a complex modulus that follows the four parameter fractional model in the frequency domain, where E r = 0.353 GPa and E u = 3.462 GPa are the relaxed and unrelaxed moduli, respectively, τ = 314.9 µs is the relaxation time and α = 0.873 is the fractional parameter. These values are taken from the characterisation of the AISI T 316L stainless steel laminated sheet and the Soundown Vibration Damping Tile material performed in [31]. A value of 0.3 is considered for the Poisson's ratio of the steel and viscoelastic layers. Three cases are studied for both FLD and CLD configurations, varying the thickness of the viscoelastic layer. For the FLD, the steel layer has a thickness H 1 = 2 mm and the thickness of the viscoelastic layer H 2 takes the value of 2 mm, 6 mm and 10 mm; for the CLD, the steel layers have a thickness H 1 = H 3 = 1 mm while the thickness of the viscoelastic layer H 2 can be 1 mm, 5 mm or 10 mm. This ensures that the homogenised model is tested for both thin and thick structures.
The process of testing the performance of the homogenised formulation goes as follows: • For each case-FLD beam, FLD plate, CLD beam and CLD plate-two finite element models are created: a beam or plate one, depending on the case, and a reference 3D model.

•
The standard and the homogenised formulations are implemented in the beam or plate model.

•
The natural frequencies, loss factor and frequency response of the structure are computed using the two finite elements models and, thus, the three formulations. For each case, three values of the thickness of the viscoelastic layer are considered.

•
The accuracy of the homogenised formulation is measured as its ability to provide results close to those obtained by a reference 3D finite element model.

Unconstrained Layer Damping (FLD)
First of all, the homogenised formulation is applied to unconstrained layer damping (FLD) beams and plates.

FLD Beams
To check the performance of the homogenised formulation for thin and thick viscoelastic layered beams, first a simply supported FLD beam 0.12 m long is considered, so that the thickness to length ratio is about 3%, 5% and 10% when the viscoelastic layer has a thickness of 2 mm, 6 mm and 10 mm, respectively.
Two models are created: a 1D model based on beam elements with 2 DOF in each end (the transverse displacement v and the rotation θ z , see Appendix A for the details) and a 3D model based on quadratic hexaedra with 27 nodes and 3 DOF in each node (the displacements u, v and w in the directions x, y and z, respectively, see [32] for the details). In the 1D model, both the Oberst and the homogenised formulation are implemented, while the 3D model serves as the reference.
In order to ensure convergence for the first 10 modes, the 1D beam model is meshed with 60 elements in its length and the 3D quadratic beam model has four elements in its width, 48 in its length and 2 in its height for each layer. This way, the beam model has 122 DOF while the solid model has 23 571 DOF, roughly 200 times more, which has a direct impact on the computation time.
The simply supported boundary conditions are modelled differently depending on the model. For the 1D model, the transverse displacement of the two ends of the beam is blocked, but the rotation is allowed. The 3D model requires additional considerations: the boundary conditions should be applied in the neutral fibre to allow the rotation of the section and the displacement in the direction of the beam should be allowed in one of the supports (Figure 3). As the position of the neutral fibre changes in frequency because damping makes it a complex magnitude, in this work, the decision to apply the boundary conditions in the middle of the section is taken. Given the stiffness of the beam, this simplification does not introduce a considerable error.  Table 1 shows the first three natural frequencies and loss factors for the three models and the three thickness cases analysed. Only the bending modes perpendicular to the beam are considered in the comparison, which are actually the ones that produce noise. In order to identify these mode shapes, the Modal Assurance Criterion (MAC) value is used, which for mode shapes φ i and φ j , takes the form where (•) H stands for the Hermitian transpose. The natural frequencies given by the homogenised formulation lie between the ones produced by the solids and those of the standard model. This means that the homogenised formulation is stiffer than the solid model because it does not introduce the in-plane displacement, but softer than the standard because the effect of shear is considered.
As expected, the Oberst method is accurate for the thinnest case producing results that diverge less than a 1% from the reference model but starts to fall short when the thickness of the viscoelastic layer increases, especially for the third mode. The homogenised method, on the contrary, is able to accurately compute the natural frequencies of the FLD beam even for the greatest thickness, with a maximum divergence of less than 4% from the values obtained with the solid model while the Oberst formulation presents a deviation of nearly 30%.
Regarding the loss factor, the values of the loss factor are not directly comparable because of the differences in natural frequencies between the three models [22]. In fact, in this case, the differences between the three models are minimal for the thinnest case but start growing together with thickness and mode number or, in other words, when the divergences in frequency increase.
On the subject of mode shapes, the standard and homogenised models produce, due to their formulation, the same eingenvectors and, thus, the MAC value for both models must be identical. In addition, both beam models are able to correctly reproduce the mode shapes of the reference model as the MAC value for every case is one.
As for the response, it is computed between 0 and 10 kHz and considering 500 samples in the frequency range so that the resolution is acceptable. Figure 4 shows the amplitude of the transverse displacement in frequency given by the three models for the three thickness cases studied. Due to the distributed pressure applied to the beam, only the symmetric modes appear in the response. In view of the results, the main conclusion is that the homogenised model is able to reproduce the response obtained by the solid model in terms of amplitude and position of the resonant peaks for every case, while the Oberst formulation misses the highest frequency peak in the thick case.

FLD Plates
The homogenised formulation is then tested in a 0.1 m × 0.1 m simply supported FLD plate. With this aim, a 2D model based on rectangular plate elements with 3 DOF in each node (the vertical displacement v and the rotations around the x and z axes θ x and θ z (see Appendix A for the details) and a 3D model based in quadratic hexaedra are implemented. Both the Oberst and homogenised formulations are implemented in the 2D model.
The 3D model has 20 elements in each direction in the plane and two in the height in each layer, while the plate model has 50 elements in each direction in the plane. This mesh ensures the convergence of the first 10 modes for the three cases. In this case, the size of the problem is eight times bigger for the solid model, with 65,559 DOF, than for the plate model that only has 7803 DOF.
The boundary conditions are applied in the four edges of the structure following the directions given in Section 4.1.1: in the plate model, transverse displacement is blocked in the edges while rotations are allowed; in the solid model, boundary conditions are applied in the middle line of the cross section and in such a manner that the displacement in the plane of the plate is allowed.
The natural frequencies and loss factors of the FLD plate for the three different thicknesses are gathered in Table 2. The second mode is double because of the symmetry of the structure. In addition, in this case, the differences in accuracy between the Oberst and the homogenised model when computing natural frequencies increase with the thickness of the viscoelastic layer, ranging from about 1% for both models when H 2 = 2 mm to 43% for the Oberst model and 7% for the homogenised when H 2 = 10 mm and the highest natural frequency. This fact is a consequence of the inclusion of the shear in the homogenised model, whose effect is more noticeable when the plate is thick.
On the subject of mode shapes, both the Oberst and the homogenised model are able to reproduce properly the results obtained in the reference model, taking into account the MAC values lie between 0.8 and 0.95, the lower value being that of the mode of highest natural frequency. This difference can be attributed to the discretisation. MAC value is omitted in Table 2 for the second mode because, being a double mode, it is not applicable.
Regarding the response in the frequency domain, similar trends as the ones found when analysing FLD beams are observed ( Figure 5): the homogenised formulation reproduces correctly the results given by the reference model while the Oberst model fails to predict the high frequency peak of the thick plate. In view of the above, it can be concluded that the homogenised model offers a clear advantage when computing the natural frequencies, mode shapes and frequency response of FLD plates, taking into account that the obtained values are closer to the reference for the three thicknesses considered. In addition, it is even more cost-effective than the 3D solid model as it is capable of obtaining similar results at a reduced computation time because of its more limited size.

Constrained Layer Damping (CLD)
The performance of the general homogenised formulation is then tested with constrained layer damping (CLD) beams and plates.

CLD Beams
The structure under study in this case is a L = 0.12 m long simply supported CLD beam. As previously stated, three cases are analysed, varying the thickness of the viscoelastic layer (1 mm, 5 mm and 10 mm) while that of the metallic layers is kept constant. This allows for testing the homogenised formulation for both thin and thick beams, as the thickness to length relationship ranges from 2.5%, for the thin case, to 10%, for the thick one.
A beam and a solid model are also developed for this case, both considering the same amount of nodes as the ones used to model FLD beam in Section 4.1.1. The single difference is that, for CLD structures, the standard formulation is the RKU that is included in Appendix B for the sake of completeness.
Boundary conditions are applied analogously as described when dealing with FLD plates: in the beam model, the transverse displacement is blocked in the two ends while the rotation is allowed; in the solid model, all displacements are blocked in the middle of the section in one end, while the displacement in the direction of the beam is allowed in the other.
The natural frequencies and loss factors computed by the three methods are shown in Table 3. In this case, the deviations in the natural frequencies are similar for the two beam models and range from less than 1% for the thin beam to about 12% for the thick beam and the highest natural frequency. Even if the differences are not as high as in the previous cases, it is worth noting that the RKU model overestimates the natural frequencies because it overestimates the equivalent stiffness [33]. Table 3. Natural frequencies and loss factors of the simply supported CLD beam computed with the three models. For the beam models, the percent difference in the natural frequencies and MAC value with respect to the 3D model are also computed. As for the mode shapes, the MAC value agrees for the RKU and homogenised formulations for all the thickness cases studied and its value decreases when the frequency of the mode increases. Apart from errors introduced by the discretisation, this is related to the coupling with bending in the other plane or the in-plane displacements that are considered in the beam model and are more relevant when either the natural frequency or the thickness of the beam are higher.
As a final consideration regarding modes, it should be taken into account that, if the viscoelastic layer is thick, the bending modes in x − z plane and the torsion modes appear at low frequency. As a beam model is not able to capture that behaviour, a plate model should be used instead.
Regarding the response (Figure 6), also in this case the divergences between the standard model and the homogenised formulation are more noticeable when either the frequency or the thickness of the viscoelastic layer increase. There is one aspect that is specific to CLD structures: the deviation of the natural frequencies computed by the homogenised formulation from those obtained by the reference 3D model can be attributed to the coupling between in-plane and out-of-plane deformation that is not included in the homogenised model.
In conclusion, the homogenised formulation outperforms the RKU for the three thickness cases under consideration, although the differences in performance are not as striking as in FLD beams because RKU was actually developed for simply supported CLD beams.

CLD Plates
The last structure analysed is a 0.1 m × 0.1 m CLD simply supported plate. Three cases are studied varying the thickness of the viscoelastic layer (1 mm, 5 mm and 10 mm) while the thickness of the metallic layers remains constant. The meshes and boundary conditions of the finite element beam and solid models are identical to those described for FLD plates in Section 4.1.2.
The natural frequencies and loss factor for the simply supported CLD plate and the three values of thickness are shown in Table 4. In addition, in this case, the second mode is double due to the symmetry of the plate. This is the case in which the homogenised model shows the greatest deviations in comparison to the reference, reaching a maximum error of 13% for the third mode of the thick plate. Even in this case, it is more accurate than the RKU that deviates about 35% from the reference model in the same conditions.
The differences between the two plate models are not so high if the mode shapes are studied, MAC values ranging around 0.73-0.95 for the cases in which it is applicable. In addition, in this case, the MAC value decreases when the frequency increases for the reason described above. Concerning the response (Figure 7), the trends are similar to those found for the CLD beam: divergences in natural frequencies increase with thickness and frequency, being especially severe for the thick plate. It should be highlighted that the RKU is not able to reproduce the high frequency peak and presents accuracy problems for the second one even for the thin plate. The homogenised formulation, on the contrary, correctly reproduces the frequency response for the thin plate in terms of amplitude and frequency. When both the thickness and the frequency increase, the deviation starts to be noticeable, even for the homogenised formulation, especially regarding the amplitude. This is related to the contribution of coupled in-plane and out-of-plane motion that is higher for CLD plates than for any of the structures analysed in this work.

Conclusions
The work presents a general homogenised formulation to study the dynamic behaviour of viscoelastic layered beam and plates by defining a frequency dependent flexural modulus that takes into account the effect of shear.
This approach is more accurate than the Oberst model for FLD structures and than the RKU formulation for CLD ones, regarding natural frequencies and frequency response, and almost reaches the accuracy of a 3D model, but at a much lower computational effort. Loss factor cannot be directly compared due to the differences in frequency.
The homogenised formulation also presents other advantages: it is valid to analyse the dynamic behaviour of viscoelastic layered structures with any number of layers; complex damping models, such as the fractional one used throughout the work, can be introduced with ease; and the modelling process is simplified because the problems derived from the locality of the boundary conditions or the variable location of the neutral fibre are avoided. On the negative side, due to its formulation, it cannot be used in the static regime and can only provide the transverse response of the structure. Given that the field of application of layered viscoelastic structures is vibration and noise reduction, these limitations should not reduce the convenience of the presented formulation.
In addition, as the presented approach allows for producing the mass and stiffness matrices of the system in an external program using conventional finite elements and include the frequency dependent stiffness modulus when computing the eigenvalues, eigenvectors or response, it can be applied to viscoelastic layered structures of any shape.
Finally, it should be noted that the homogenised formulation presented can be applied to both thin and thick viscoelastic layered structures, but it is specially advantageous in this last case for which the only available solutions are computationally expensive.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Elementary Matrices
Appendix A. 1

. Beam Matrices
The mass matrix of the beam element used throughout the work and shown in Figure A1 is A and l are the transverse section of the beam and the length of the element, respectively. It should be noted that, even if consistent mass matrices are used in the work, the presented formulation is also applicable if the mass matrix is lumped. Using lumped mass matrices could, in fact, be advantageous in terms of computation time when solving the eigenvalue and eigenvector problem or when dealing with transitory dynamics, but at the expense of a loss in accuracy.

Appendix A.2. Plate Matrices
The consistent mass matrix of the plate element used throughout the work and shown in Figure A2 is  Figure A2. Plate element of size a × b with four nodes and 3 DOF in each node, the transverse displacement v and the rotations θ x and θ z .

Appendix B. RKU Formulation
The Ross, Kerwin and Ungar formulation defines the equivalent flexural stiffness for three layered viscoelastic beams as the flexural stiffnesses for the ith layer being

(A6)
and where H i is the thickness of the ith layer and G i its shear modulus. The variable k * B stands for the wave number of the bending wave and its defined as For plates, the flexural stiffness D * eq and the mass per unit area ρ S should be used instead and, as the computations are done for unit width, b = 1 units of length (in the considered coherent system of units).