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

Abstract

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.

1. 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 dealing with simple structures and damping models, but, when the geometry of the structure under study, its boundary conditions or applied forces are complex, analytical models fall short and finite element models are the common choice. Two paths can be followed if solid finite element models that are demanding in terms of computational resources [7] are to be avoided: to propose a homogenised formulation that takes into account the contribution of every layer so that traditional plate elements can be used; or to develop a specific finite element model for layered structures. In the engineering practise, it is also common to use the multilayer elements available in commercial finite element software, but they do not allow the introduction of complex damping models that are often needed when dealing with these kinds of surface treatments.

Models can be divided according to the technology for which they have been developed. Regarding the FLD case, homogenised formulations have been commonly proposed, starting from the Oberst model [8], in which an equivalent flexural stiffness for the two layers is proposed and up to the more recent [6], in which a model of a cantilever plate was developed with the aim of studying the frequency dependence of the response. Neither of these homogenised models considered the effect of shear and are thus limited to thin plates where this effect can be neglected.

As for CLD structures, the first works aiming to understand their dynamic behaviour are the ones by Ross [9], Kerwin [10] and Ungar [11] that resulted in the RKU method that is still nowadays the standard in the field. This method defines an equivalent flexural stiffness that links the shear strain in the damping layer to the bending motion of the plate. The analytical procedure has been followed in subsequent works such as [12], in which the work by Ross, Kerwin and Ungar was extended or, after, in [13] in which a model for thick layered structures that considers nine degrees of freedom (DOF) is presented. This last work is stated as a general case and can degenerate into one- or two-layered structures so it can also be applied to FLD structures. Due to the complexity of the analytical modelling, these models are limited to simple geometries.

In addition, in [14], a homogenised formulation similar to those for modelling FLD structures is presented. Thin plate elements that consider the properties of the three layers are proposed but, as shear is not considered, it is only valid when the viscoelastic core is thin.

With the idea of dealing with complex geometries, applied forces or boundary conditions, different finite element models have been proposed. For example, in [15], an 8 DOF beam element is presented; in [16], triangular and rectangular sandwich elements with seven DOF in each node are developed; in [17], a three-layer four-node rectangular element with 7 DOF on every node is used, or, in [18], a nine-node isoparametric 2D element is proposed to model the behaviour of active-passive composite beams. The finite element method is also used, but with a different approach in the series [19,20], where the viscoelastic core of a CLD structure is modelled as a thick beam or plate, depending on the structure under study, and the metallic layers are represented as thin beams or plates. The main drawback of these methods is the need for developing a specific finite element formulation for the application, instead of taking advantage of the solutions already available.

In view of the above, this work presents a general homogenised formulation for the dynamic analysis of thick viscoelastic layered beams and plates that extends the previous models proposed by the authors [21,22,23,24] to N-layered structures. This generalisation addresses the need for a single approach valid for both one- or two-dimensional models, regardless of the type of viscoelastic surface treatment under study and its thickness, and that can deal with the frequency dependent behaviour of the viscoelastic layer.

The formulation uses conventional beam and plate finite elements and introduces the effect of shear by means of a frequency dependent equivalent flexural stiffness. With this aim, first, the formulation for a general one- or two-dimensional surface treatment with any number of layers is presented; then, the method is implemented in a finite element model and applied to FLD and CLD beams and plates for which the eigenpairs, i.e., eigenvalues and eigenvectors, and the dynamic response to an external force are computed; finally, the obtained results are compared to the ones provided by a three-dimensional finite element model and to the standard model in the field: the Oberst model for FLD structures and 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.

2. 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_{\mathrm{eq}}$ for a N-layered beam can be expressed as the sum of the flexural stiffness of its layers as

$$B_{\mathrm{eq}}=\sum _{i=1}^N{B}_i=\sum _{i=1}^N{E}_i{I}_i,$$

(1)

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

$$I_i=b{\int }_{h_{i-1}}^{h_i}{y}^2\mathrm{d}y,$$

(2)

where b is the width of the beam that is supposed to be constant, and defining the integration limits

$$\begin{array}{cc}\hfill & h_0=-h_{\mathrm{n}}\hfill \\ \hfill & h_1=-h_{\mathrm{n}}+H_1\hfill \\ \hfill & \dots \hfill \\ \hfill & h_i=-h_{\mathrm{n}}+H_1+\dots +H_i\hfill \\ \hfill & \dots \hfill \\ \hfill & h_{\mathrm{N}}=-h_{\mathrm{n}}+H_1+\dots +H_N,\hfill \end{array}$$

(3)

for every layer, where $h_{\mathrm{n}}$ stands for the position of the neutral fibre that is located in

$$h_{\mathrm{n}}=\frac{E_1{H}_1\frac{H_1}{2}+E_2{H}_2\left(H_1+\frac{H_2}{2}\right)+\dots +E_N{H}_N\left(H_1+H_2+\dots +\frac{H_N}{2}\right)}{{\displaystyle \sum _{i=1}^N{E}_i{H}_i}}.$$

(4)

If the small deformation hypothesis holds, it is possible to uncouple the transverse displacement due to bending $v_{\mathrm{M}}(x,t)$ and to shear $v_{\mathrm{Q}}(x,t)$ as

$$v(x,t)=v_{\mathrm{M}}(x,t)+v_{\mathrm{Q}}(x,t).$$

(5)

The first term, as it is originated by the bending moment $M(x,t)$, satisfies

$$M(x,t)=B_{\mathrm{eq}}\frac{{\partial }^2{v}_{\mathrm{M}}(x,t)}{\partial x^2},$$

(6)

where $B_{\mathrm{eq}}$ is the equivalent flexural stiffness from (1).

Likewise, the second term is related to the shear force $Q(x,t)$ by

$$Q(x,t)=K_{\mathrm{eq}}\frac{\partial v_{\mathrm{Q}}(x,t)}{\partial x},$$

(7)

where the total shear stiffness $K_{\mathrm{eq}}$ is obtained by combining the shear stiffnesses of all layers

$$\frac{1}{K_{\mathrm{eq}}}=\sum _{i=1}^N\frac{1}{K_i}.$$

(8)

If a second order shear stress law is considered, the shear stiffness for the ith layer follows

$$\frac{1}{K_i}=\frac{1}{B_{\mathrm{eq}}^{2}b}{\int }_{h_{i-1}}^{h_i}\frac{{\Gamma }_i{\left(y\right)}^2}{G_i}\mathrm{d}y,$$

(9)

where $G_i$ is the shear modulus of the material of the ith layer and

$$\Gamma \left(y\right)={\int }_{S^{\prime }}yE\left(y\right)\mathrm{d}{S}^{\prime },$$

(10)

where $S^{\prime }$ 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), ${\Gamma }_i\left(y\right)$ has the value

$${\Gamma }_1\left(y\right)=\frac{E_{1}b\left(h_0^2-y^2\right)}{2}$$

(11)

for the bottom layer,

$${\Gamma }_i\left(y\right)=\sum _{j=1}^{i-1}\frac{E_{j}b\left(h_{j-1}^2-h_j^2\right)}{2}+\frac{E_{i}b\left(h_{i-1}^2-y^2\right)}{2}$$

(12)

for the middle layers 2 to N-1 and

$${\Gamma }_N\left(y\right)=\frac{E_{N}b\left(h_N^2-y^2\right)}{2}$$

(13)

for the top layer N.

The equivalent flexural and shear stiffnesses $B_{\mathrm{eq}}$ and $K_{\mathrm{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

$$\frac{B_{\mathrm{eq}}}{{\rho }_{\mathrm{L}}}\frac{{\partial }^{4}v(x,t)}{\partial x^4}-\frac{B_{\mathrm{eq}}}{K_{\mathrm{eq}}}\frac{{\partial }^{4}v(x,t)}{\partial x^2\partial t^2}+\frac{{\partial }^{2}v(x,t)}{\partial t^2}=q(x,t),$$

(14)

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 ${\rho }_{\mathrm{L}}={\sum }_{i=1}^N{\rho }_i{A}_i$ being the mass per unit length and ${\rho }_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 the form

$$v(x,t)=A\exp \left[\mathrm{i}(\beta x-\omega t)\right],$$

(15)

where A, $\beta$ and $\omega$ 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

$$\frac{\tilde{B}\left(t\right)}{{\rho }_{\mathrm{L}}}\frac{{\partial }^{4}v(x,t)}{\partial x^4}+\frac{{\partial }^{2}v(x,t)}{\partial t^2}=q(x,t),$$

(16)

where $\tilde{B}\left(t\right)$ is the time domain counterpart of the flexural stiffness modulus that, in the frequency domain, takes the form

$$B\left(\omega \right)=\frac{B_{\mathrm{eq}}}{{\left(\phi \left(\omega \right)+\sqrt{\phi {\left(\omega \right)}^2+1}\right)}^2},$$

(17)

in which

$$\phi \left(\omega \right)=\frac{\omega \sqrt{B_{\mathrm{eq}}{\rho }_{\mathrm{L}}}}{2K_{\mathrm{eq}}}.$$

(18)

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 ${\rho }_{\mathrm{L}}$ should be substituted by the mass per unit area

$${\rho }_{\mathrm{S}}=\sum _{i=1}^N{\rho }_i{H}_i,$$

(19)

where ${\rho }_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_{\mathrm{eq}}$ should be substituted by its counterpart for plates

$$D_{\mathrm{eq}}=\sum _{i=1}^N{\int }_{h_{i-1}}^{h_i}\frac{E_i}{3(1-{\nu }_i^2)}{z}^2\mathrm{d}z,$$

(20)

$E_i$ and ${\nu }_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

$$D\left(\omega \right)=\frac{D_{\mathrm{eq}}}{{\left(\phi \left(\omega \right)+\sqrt{\phi {\left(\omega \right)}^2+1}\right)}^2},$$

(21)

where $\phi \left(\omega \right)$ follows the expression

$$\phi \left(\omega \right)=\frac{\omega \sqrt{D_{\mathrm{eq}}{\rho }_{\mathrm{S}}}}{2K_{\mathrm{eq}}}.$$

(22)

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 $\omega$, the Kirchhoff–Love thin plate equation [27] becomes

$$\tilde{D}\left(t\right)\left(\frac{{\partial }^{4}v(x,z,t)}{\partial x^4}+2\frac{{\partial }^{4}v(x,z,t)}{\partial x^2\partial z^2}+\frac{{\partial }^{4}v(x,z,t)}{\partial z^4}\right)+{\rho }_{\mathrm{S}}\frac{{\partial }^{2}v(x,z,t)}{\partial t^2}=p(x,z,t),$$

(23)

where $p(x,z,t)$ is the distributed force per unit area applied to the plate and $\tilde{D}\left(t\right)$ 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.

3. 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 $\mathbf{K}\left(0\right)$ as

$${\mathbf{K}}_{\mathrm{beam}}^{*}\left(\omega \right)=\frac{B^{*}\left(\omega \right)}{B\left(0\right)}\mathbf{K}\left(0\right)$$

(24)

for beams and

$${\mathbf{K}}_{\mathrm{plate}}^{*}\left(\omega \right)=\frac{D^{*}\left(\omega \right)}{D\left(0\right)}\mathbf{K}\left(0\right)$$

(25)

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.

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

$$\left(-{\lambda }_r^{*}\mathbf{M}+{\mathbf{K}}^{*}\left({\omega }_r\right)\right){\mathsf{\varphi }}_r^{*}=\mathbf{0},$$

(26)

where $\mathbf{M}$ is the mass matrix, ${\omega }_r=\mathrm{Re}\left(\sqrt{{\lambda }_r^{*}}\right)$ is the natural frequency for mode rth and ${\lambda }_r^{*}$ and ${\mathsf{\varphi }}_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 $\mathbf{K}\left(0\right)$ 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

$${\eta }_r=\frac{\mathrm{Im}\left({\lambda }_r^{*}\right)}{\mathrm{Re}\left({\lambda }_r^{*}\right)}.$$

(27)

3.2. Dynamic Response

The dynamic response of the structure is obtained by solving the motion equation in the frequency domain

$$\left(-{\omega }^2\mathbf{M}+{\mathbf{K}}^{*}\left(\omega \right)\right){\mathbf{v}}^{*}\left(\omega \right)=\mathbf{F}\left(\omega \right)$$

(28)

for the response vector ${\mathbf{v}}^{*}\left(\omega \right)$.

For the computations, a uniform impulse pressure (1 Pa) on the upper surface of the structure is considered as the force vector $\mathbf{F}\left(\omega \right)$ and, for the comparison between models, the RMS value of the transverse displacement ${\mathbf{v}}^{*}\left(\omega \right)$ on that same surface is computed as

$$v_{\mathrm{RMS}}\left(\omega \right)=\sqrt{\frac{1}{R}\sum _{i=1}^R{\left|v_i^{*}\left(\omega \right)\right|}^2},$$

(29)

where R is the number of nodes of the model and $v_i^{*}\left(\omega \right)$ the response in the frequency domain of the ith node in the transverse direction.

4. 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 ${\rho }_1$ = 7782 kg/m³ and their Young’s modulus $E_1$ = 176.24 GPa. The viscoelastic layer has a density of ${\rho }_2$ = 1423 kg/m³ and a complex modulus that follows the four parameter fractional model

$$E_2^{*}\left(\omega \right)=\frac{E_{\mathrm{r}}+E_{\mathrm{u}}{\left(\mathrm{i}\omega \tau \right)}^{\alpha }}{1+{\left(\mathrm{i}\omega \tau \right)}^{\alpha }}$$

(30)

in the frequency domain, where $E_{\mathrm{r}}$ = 0.353 GPa and $E_{\mathrm{u}}$ = 3.462 GPa are the relaxed and unrelaxed moduli, respectively, $\tau$ = 314.9 $\mathsf{\mu }$s is the relaxation time and $\alpha$ = 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.

4.1. Unconstrained Layer Damping (FLD)

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

4.1.1. 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 ${\theta }_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. Natural frequencies and loss factors of the simply supported FLD beam computed with the three models. For the beam models, the  percent difference in the natural frequencies and MAC value in comparison to the 3D model are shown.

${\mathit{H}}_2$ (mm)	Mode	3D Model		Homogenised Model				Oberst Model
		$\mathit{\omega }$ (rad/s)	$\mathit{\eta }$	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC
2	1	1813	0.0760	1813	0.001	0.0765	1.00	1814	0.05	0.0766	1.00
	2	7578	0.0600	7578	0.006	0.0602	1.00	7596	0.24	0.0599	1.00
	3	17,201	0.0345	17,204	0.02	0.0348	1.00	17,299	0.57	0.0346	1.00
6	1	2450	0.4472	2453	0.25	0.4520	1.00	2482	1.32	0.4502	1.00
	2	11,371	0.1839	11,467	0.87	0.1904	1.00	11,917	4.80	0.1819	1.00
	3	25,035	0.0953	25,369	1.33	0.1032	1.00	27,501	9.85	0.0930	1.00
10	1	4013	0.4778	4056	1.08	0.4916	1.00	4208	4.86	0.4663	1.00
	2	16,683	0.1655	17,128	2.67	0.1818	1.00	19,204	15.1	0.1491	1.00
	3	33,861	0.0895	35,093	3.64	0.1096	1.00	43,905	29.7	0.0745	1.00

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 ${\mathsf{\varphi }}_i$ and ${\mathsf{\varphi }}_j$, takes the form

$$\mathrm{MAC}=\frac{|{\mathsf{\varphi }}_i^{\mathrm{H}}{\mathsf{\varphi }}_j{|}^2}{|{\mathsf{\varphi }}_i^{\mathrm{H}}{\mathsf{\varphi }}_i\left|\right|{\mathsf{\varphi }}_j^{\mathrm{H}}{\mathsf{\varphi }}_j|},$$

(31)

where ${\left(•\right)}^{\mathrm{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.

4.1.2. 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 ${\theta }_x$ and ${\theta }_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. Natural frequencies and loss factors of the simply supported FLD plate computed with the three models. For the plate models, the  percent difference in frequency and MAC value in comparison to the 3D models is shown. MAC value is omitted for double modes.

${\mathit{H}}_2$ (mm)	Mode	3D Model		Homogenised Model				Oberst Model
		$\mathit{\omega }$ (rad/s)	$\mathit{\eta }$	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC
2	1	5631	0.0683	5681	0.88	0.0692	0.95	5691	1.07	0.0695	0.95
	2	14,318	0.0392	14,408	0.68	0.0397	-	14,477	1.11	0.0401	-
	3	22,837	0.0272	23,064	0.99	0.0277	0.80	23,242	1.77	0.0281	0.80
6	1	8225	0.2326	8565	4.14	0.2353	0.95	8836	7.43	0.2273	0.95
	2	20,722	0.1110	21,383	3.17	0.1182	-	22,976	10.9	0.1079	-
	3	31,868	0.0754	33,182	4.12	0.0837	0.80	37,027	16.2	0.0725	0.80
10	1	12,072	0.2171	13,048	8.09	0.2223	0.95	14,334	18.7	0.1893	0.95
	2	28,186	0.1051	30,001	6.44	0.1299	-	36,725	30.3	0.0867	-
	3	41,125	0.0746	44,099	7.23	0.0949	0.80	59,023	43.5	0.0578	0.80

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

4.2. Constrained Layer Damping (CLD)

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

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

${\mathit{H}}_2$ (mm)	Mode	3D Model		Homogenised Model				RKU Model
		$\mathit{\omega }$ (rad/s)	$\mathit{\eta }$	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC
1	1	3162	0.0337	3165	0.12	0.0345	1.00	3166	0.14	0.0340	1.00
	2	11,934	0.0379	11,959	0.21	0.0421	1.00	11,984	0.42	0.0407	1.00
	3	24,741	0.0370	24,752	0.05	0.0459	1.00	24,916	0.71	0.0429	1.00
5	1	7185	0.0822	7301	1.62	0.0869	1.00	7283	1.36	0.0832	1.00
	2	22,746	0.0733	23,448	3.09	0.0953	1.00	23,474	3.20	0.0921	1.00
	3	40,341	0.0582	41,807	3.63	0.0928	0.99	42,105	4.37	0.0882	0.99
10	1	10,131	0.1096	10,707	5.69	0.1166	1.00	10,640	5.02	0.1099	1.00
	2	27,963	0.0846	30,743	9.95	0.1162	0.97	30,728	9.89	0.1124	0.97
	3	46,179	0.0654	51,473	11.5	0.1035	0.93	51,723	12.0	0.0998	0.93

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

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.

4.2.2. 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. Natural frequencies and loss factors of the simply supported CLD plate computed with the three models. For the plate models, the  percent difference in natural frequency and MAC values in comparison to the 3D model are shown. The MAC value is omitted for double modes.

${\mathit{H}}_2$ (mm)	Mode	3D Model		Homogenised Model				RKU Model
		$\mathit{\omega }$ (rad/s)	$\mathit{\eta }$	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC	$\mathit{\omega }$ (rad/s)	Diff (%)	$\mathit{\eta }$	MAC
1	1	9125	0.0554	9162	0.41	0.0426	0.95	10,527	15.4	0.0337	0.95
	2	20,737	0.0420	21,103	1.76	0.0473	-	24,480	18.1	0.0364	-
	3	30,906	0.0393	31,459	1.36	0.0492	0.80	36,819	19.1	0.0368	0.80
5	1	17,786	0.0877	18,633	4.77	0.0987	0.95	21,692	22.0	0.0782	0.95
	2	34,684	0.0667	36,547	5.37	0.0977	-	43,215	24.6	0.0773	-
	3	47,244	0.0675	49,696	5.09	0.0929	0.78	59,319	25.6	0.0728	0.78
10	1	22,901	0.1040	24,930	12.0	0.1229	0.94	29,216	31.2	0.0979	0.94
	2	40,116	0.0739	45,397	13.5	0.1106	-	53,953	34.9	0.0893	-
	3	53,087	0.0608	59,928	13.0	0.0999	0.73	71,706	35.2	0.0805	0.73

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

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