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

Authors: Fernando Cortés; Imanol Sarría
Journal: Shock and Vibration
 (2015)

---

## Abstract

This paper presents an analysis of the dynamic behaviour of constrained layer damping (CLD) beams with thick viscoelastic layer. A homogenised model for the flexural stiffness is formulated using Reddy-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 2D model 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.


---

## Body

## 1. 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 Figure1 the cross section of a three-layer composed beam is represented, the thickness of the base, viscoelastic, and constraining layers being H
1, H
2, and H
3, respectively, and b being the width of the beam. Although in CLD beams H
2
≪
H
1, in this work three-layer sandwich beams with any value for H
2 are considered. The properties for these elastic materials are Young’s moduli E
1 and E
3, Poisson’s coefficients ν
1 and ν
3, and densities ρ
1 and ρ
3, whereas those of the viscoelastic core are E
2
∗, ν
2, and ρ
2, respectively. If the materials are considered isotropic, the shear moduli G
1, G
2
∗, and G
3 are related to the extensional ones by G
1
=
E
1
/
2
(
1
+
ν
1
), G
2
∗
=
E
2
∗
/
2
(
1
+
ν
2
), and G
3
=
E
3
/
2
(
1
+
ν
3
) for the three layers.Figure 1
Cross section of a three-layer sandwich beam with a viscoelastic damping layer.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:(1)
η
2
=
E
2
″
E
2
′
,where E
2
′ and E
2
″ are the real and imaginary components of the complex modulus E
2
∗
=
E
2
′
+
i
E
2
″, which are known as the storage and loss modulus, respectively, the former being also frequently represented by E
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 self-supporting 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–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.


## 2. 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 stiffnessB
eq may be obtained by the addition of the individual contribution of each layer:(2)
B
eq
=
B
1
+
B
2
+
B
3
,where B
1
=
E
1
I
1, B
2
=
E
2
I
2, and B
3
=
E
3
I
3 are the complex flexural stiffness of the three layers and I
1, I
2, and I
3 are the complex cross-sectional second order moments, computed with respect to the neutral axis, given by (3)
I
1
=
1
12
b
H
1
3
+
b
H
1
h
n
-
H
1
2
2
,
(4)
I
2
=
1
12
b
H
2
3
+
b
H
2
H
1
+
H
2
2
-
h
n
2
,
(5)
I
3
=
1
12
b
H
3
3
+
b
H
3
H
1
+
H
2
+
H
3
2
-
h
n
2
,respectively, in which the complex position of the neutral axis is represented by (6)
h
n
=
E
1
H
1
2
/
2

+

E
2
H
2
(
H
1

+

H
2
/
2
)

+

E
3
H
3
(
H
1
+
H
2
+
H
3
/
2
)
E
1
H
1
+
E
2
H
2
+
E
3
H
3
.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 (7)
1
K
eq
=
1
K
1
+
1
K
2
+
1
K
3
.For the geometry represented in Figure 1, the stiffness of K
1, K
2, and K
3 of the individual layers satisfies(8)
1
K
1
=
1
G
1
B
eq
2
b
∫
-
h
n
-
h
n
+
H
1
Ω
1
2
(
y
)
d
y
,
(9)
1
K
2
=
1
G
2
B
eq
2
b
∫
-
h
n
+
H
1
-
h
n
+
H
1
+
H
2
Ω
2
2
(
y
)
d
y
,
(10)
1
K
3
=
1
G
3
B
eq
2
b
∫
-
h
n
+
H
1
+
H
2
-
h
n
+
H
1
+
H
2
+
H
3
Ω
3
2
(
y
)
d
y
,respectively, where(11)
Ω
1
2
(
y
)
=
b
2
E
1
2
h
n
2
-
y
2
2
4
,
(12)
Ω
2
2
(
y
)
=
b
2
E
1
H
1
h
n
-
H
1
2
+
E
2
h
n
-
H
1
2
-
y
2
2
2
,
(13)
Ω
3
2
(
y
)
=
b
2
E
3
2
H
1
+
H
2
+
H
3
-
h
n
2
-
y
2
2
4
,respectively. In these equations, B
eq is the flexural stiffness given by (2).By solving these integrals, it yields(14)
1
K
1
=
6
5
G
1
b
H
1
10
r
n
2
-
15
r
n
+
6
1
+
r
2
M
2
+
r
3
M
3
2
1
+
3
r
n
-
1
2
2
,
(15)
1
K
2
=
1
G
2
b
H
1
·
36
T
2
r
n
-
1
2
-
12
M
2
T
2
2
r
n
-
1
(
2
T
2
-
3
r
n
+
6
)
1
+
r
2
M
2
+
r
3
M
3
2
1
+
3
r
n
-
1
2
2
+
6
5
G
2
b
H
2
·
M
2
2
T
2
4
10
r
n
2
-
15
r
n
T
2
-
40
r
n
+
6
T
2
2
+
30
T
2
+
40
1
+
r
2
M
2
+
r
3
M
3
2
1
+
3
r
n
-
1
2
2
,
(16)
1
K
3
=
6
M
3
2
T
3
4
5
G
3
b
H
3
×
10
r
n
2
-
25
r
n
T
3
-
40
r
n
T
2
-
40
r
n
+
16
T
3
2
+
50
T
3
T
2
+
50
T
3
+
40
T
2
2
+
80
T
2
+
40
·
1
+
r
2
M
2
+
r
3
M
3
2
1
+
3
r
n
-
1
2
2
-
1
,respectively, where M
2
=
E
2
/
E
1, M
3
=
E
3
/
E
1, T
2
=
H
2
/
H
1, T
3
=
H
3
/
H
1
r
n
=
2
h
n
/
H
1, r
2
=
I
2
/
I
1, and r
3
=
I
3
/
I
1. For base and constraining layers made of metallic materials, which are much thinner and more rigid than the viscoelastic layer, the terms 1
/
K
1 and 1
/
K
3 may be neglected with respect to the coefficient of the polymeric layer 1
/
K
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:(17)
1
K
eq
=
6
5
G
S
,where S is the total cross-sectional area.By considering the shear coefficientK
eq, the flexural field equation in free vibration is given by Timoshenko’s formula [24]:(18)
B
eq
ρ
l
∂
4
v
(
x
,
t
)
∂
x
4
-
B
eq
K
eq
∂
4
v
(
x
,
t
)
∂
x
2
∂
t
2
+
∂
2
v
(
x
,
t
)
∂
t
2
=
0
,where ρ
l 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 B
K, considering shear effects and satisfying the homogenised Euler-Bernoulli field equation, can be found to be(19)
B
K
=
B
eq
1
+
φ
2
(
ω
)
+
φ
(
ω
)
2
.The function φ
(
ω
) takes into account the shear effects and is given by(20)
φ
(
ω
)
=
ω
ρ
l
B
eq
2
K
eq
.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 B
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 stiffnessK
eq tends to infinity, that is, when shear deformations are negligible. In this case, the function φ
(
ω
) tends to zero and the equivalent flexural stiffness B
K tends to the classic B
eq.


## 3. 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 isl
=
200
mm, the width is b
=
20
mm, the thickness of the base and constraining metallic layers is H
1
=
H
3
=
1
mm, and for the viscoelastic layer H
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 E
1
=
176.2
×
1
0
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 four-parameter fractional model [27, 28], given by(21)
E
2
∗
(
ω
)
=
E
r
+
E
u
i
τ
ω
α
1
+
i
τ
ω
α
,where E
r and E
u represent the relaxed and unrelaxed modulus, respectively; τ is the relaxation time and α is the fractional parameter. The parameter values are summarised in Table 1.Table 1
Parameters of the fractional derivative model.
E
r (GPa)
E
u (GPa)
τ (1
0
-
6 s)
α
0.353
3.462
314.9
0.873The 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–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.Figure 2
Finite element model for the 2D dynamical analysis of a sandwich beam.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.Figure 3
Rectangular Finite Element (Plane-Stress Assumption).Consistent Mass Matrix
M. Consider(22)
ρ
a
b
t
36
4
0
2
0
1
0
2
0
4
0
2
0
1
0
2
4
0
2
0
1
0
4
0
2
0
1
4
0
2
0
4
0
2
4
0
4
.Stiffness Matrix
K
Obtained by Reduced Integration with Hourglass Control. Consider(23)
E
t
24
γ
1
-
ν
2
×
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
3
γ
(
1
+
ν
)
3
(
1
-
ν
)
-
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
-
3
ν
)
-
3
(
1
-
ν
)
-
2
γ
2
(
2
+
ν
2
)
-
3
γ
(
1
+
ν
)
-
3
(
1
-
ν
)
+
2
γ
2
(
2
+
ν
2
)
3
γ
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)
3
γ
(
1
-
3
ν
)
-
3
γ
2
(
1
-
ν
)
+
2
(
2
+
ν
2
)
-
3
γ
(
1
+
ν
)
-
3
γ
2
(
1
-
ν
)
-
2
(
2
+
ν
2
)
-
3
b
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
-
2
(
4
-
ν
2
)
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
+
ν
)
-
3
(
1
-
ν
)
+
2
γ
2
(
2
+
ν
2
)
-
3
γ
(
1
-
3
ν
)
-
3
(
1
-
ν
)
-
2
γ
2
(
2
+
ν
2
)
3
γ
(
1
+
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)
3
γ
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
-
2
(
4
-
ν
2
)
3
γ
(
1
+
ν
)
-
3
γ
2
(
1
-
ν
)
-
2
(
2
+
ν
2
)
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
3
γ
(
1
+
ν
)
3
(
1
-
ν
)
-
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)
3
γ
(
1
-
3
ν
)
-
3
γ
2
(
1
-
ν
)
+
2
(
2
+
ν
2
)
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
+
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)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 whoseith mass M
i and stiffness K
i matrices are(24)
M
i
=
ρ
1
H
1
+
ρ
2
H
2
+
ρ
3
H
3
b
a
i
420
·
156
22
a
i
54
-
13
a
i
22
a
i
4
a
i
2
13
a
i
-
3
a
i
2
54
13
a
i
156
-
22
a
i
-
13
a
i
-
3
a
i
2
-
22
a
i
4
a
i
2
,
(25)
K
i
=
B
a
i
3
12
6
a
i
-
12
6
a
i
6
a
i
4
a
i
2
-
6
a
i
2
a
i
2
-
12
-
6
a
i
12
-
6
a
i
6
a
i
2
a
i
2
-
6
a
i
4
a
i
2
,respectively, where a
i is the length of the ith finite element. The complex flexural stiffness B of (25) is given by (2) B
=
B
eq
∗ for the RKU model and by (19) B
=
B
K
∗ for the new thick beam model. The discretisation is also made with 60 finite elements along span.


### 3.2. Extraction of Eigenvalues

The equation from which the complex eigenvalues of the system under study can be obtained is(26)
-
λ
r
∗
M
+
K
∗
ω
r
ϕ
r
∗
=
0
,where λ
r
∗ and ϕ
r
∗ are the complex eigenvalue and eigenvector of the rth mode, respectively, M is the mass matrix, and K
∗ is the complex stiffness matrix, which is dependent on frequency. This ω
r is the real part of the square root of the complex eigenvalue λ
r
∗:(27)
ω
r
=
Re
λ
r
∗
,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(28)
-
λ
r
,
0
M
+
K
∗
0
ϕ
r
,
0
=
0
,and then the undamped eigensolutions λ
r
,
0 and ϕ
r
,
0 can be obtained. Then, Nelson’s method [35] is used in order to calculate the eigenvector derivative ϕ
r
′. From this eigenvector derivative and taking into account the variation of the complex stiffness for the obtained eigenfrequency,(29)
Δ
K
∗
ω
r
,
0
=
K
∗
ω
r
,
0
-
K
∗
(
0
)
,and by means of Taylor’s series approach, a complex finite increment Δ
ϕ
r
∗ of the eigenvector is obtained. The complex eigenvector can be approximated with(30)
ϕ
r
∗
=
ϕ
r
,
0
+
Δ
ϕ
r
∗
,with which the complex eigenvalue λ
r
∗ is estimated according to(31)
λ
r
∗
=
ϕ
r
∗
H
Κ
∗
(
ω
r
)
ϕ
r
∗
ϕ
r
∗
H
M
ϕ
r
∗
,where (
·
)
H denotes the Hermitian transpose operator, that is, the complex conjugate transposition. Equations (29)–(31) can be iterated making use of the new eigenfrequency ω
r 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ω
-
r derived from (27) and loss factor derived from(32)
λ
r
∗
=
ω
-
r
2
1
+
i
η
r
.can be computed. The corresponding results for the three thicknesses and for the three models under study are shown in Tables 2–4.Table 2
Modal properties of the sandwich cantilever beam withH
2
=
1 mm.
ω
-
1 (rad/s)
η
1
ω
-
2 (rad/s)
η
2
ω
-
3 (rad/s)
η
3
2D model
411.1
0.024
2523
0.049
6851
0.047
Present model
412.2
0.019
2547
0.037
6912
0.043
RKU model
412.4
0.018
2552
0.034
6952
0.037Table 3
Modal properties of the sandwich cantilever beam withH
2
=
5 mm.
ω
-
1 (rad/s)
η
1
ω
-
2 (rad/s)
η
2
ω
-
3 (rad/s)
η
3
2D model
1003
0.080
5659
0.113
13987
0.091
Present model
1013
0.066
5909
0.094
14504
0.103
RKU model
1019
0.052
6051
0.068
15461
0.063Table 4
Modal properties of the sandwich cantilever beam withH
2
=
10 mm.
ω
-
1 (rad/s)
η
1
ω
-
2 (rad/s)
η
2
ω
-
3 (rad/s)
η
3
2D model
1561
0.118
8106
0.138
18828
0.099
Present model
1589
0.098
8714
0.127
19779
0.131
RKU model
1613
0.069
9291
0.076
22875
0.064It can be pointed out that the results for the natural frequencyω
r 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 Table3), 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η
r, 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 η
r depends on frequency according to (21), and the natural frequencies for the models are not the same. Instead of modal loss factor η
r, the amplitudes of the resonance peaks will be compared in the next section.


### 3.3. Analysis of the Frequency Response Function

The matrix system for a cantilever beam is given by(33)
-
ω
k
2
M
+
K
∗
ω
k
U
k
∗
=
F
∗
,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
k
∗ at each desired frequency ω
k. 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 F
=
1
×
exp
⁡
(
i
ω
t
)
N is applied on the right-hand side of our system.Figure 4
Frequency response function for the free-end for the cantilever beam: (a)H
2
=
1
mm, (b) H
2
=
5
mm, and (c) H
2
=
10
mm.
(a)
(b)
(c)Figure4(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, Tables5–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.Table 5
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thicknessH
2
=
1 mm.
A
1
A
2
A
3
2D model
-
6.357
-
8.608
-
10.58
Present model
-
6.299
-
8.304
-
10.46
RKU model
-
6.203
-
8.211
-
10.31Table 6
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thicknessH
2
=
5 mm.
A
1
A
2
A
3
2D model
-
7.552
-
11.44
-
13.10
Present model
-
7.285
-
11.16
-
13.04
RKU model
-
7.058
-
10.87
-
12.68Table 7
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thicknessH
2
=
10 mm.
A
1
A
2
A
3
2D model
-
9.120
-
12.71
-
14.15
Present model
-
8.812
-
12.48
-
14.13
RKU model
-
8.451
-
12.06
-
13.72Since 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.


## 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 isl
=
200
mm, the width is b
=
20
mm, the thickness of the base and constraining metallic layers is H
1
=
H
3
=
1
mm, and for the viscoelastic layer H
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 E
1
=
176.2
×
1
0
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 four-parameter fractional model [27, 28], given by(21)
E
2
∗
(
ω
)
=
E
r
+
E
u
i
τ
ω
α
1
+
i
τ
ω
α
,where E
r and E
u represent the relaxed and unrelaxed modulus, respectively; τ is the relaxation time and α is the fractional parameter. The parameter values are summarised in Table 1.Table 1
Parameters of the fractional derivative model.
E
r (GPa)
E
u (GPa)
τ (1
0
-
6 s)
α
0.353
3.462
314.9
0.873The 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–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.Figure 2
Finite element model for the 2D dynamical analysis of a sandwich beam.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.Figure 3
Rectangular Finite Element (Plane-Stress Assumption).Consistent Mass Matrix
M. Consider(22)
ρ
a
b
t
36
4
0
2
0
1
0
2
0
4
0
2
0
1
0
2
4
0
2
0
1
0
4
0
2
0
1
4
0
2
0
4
0
2
4
0
4
.Stiffness Matrix
K
Obtained by Reduced Integration with Hourglass Control. Consider(23)
E
t
24
γ
1
-
ν
2
×
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
3
γ
(
1
+
ν
)
3
(
1
-
ν
)
-
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
-
3
ν
)
-
3
(
1
-
ν
)
-
2
γ
2
(
2
+
ν
2
)
-
3
γ
(
1
+
ν
)
-
3
(
1
-
ν
)
+
2
γ
2
(
2
+
ν
2
)
3
γ
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)
3
γ
(
1
-
3
ν
)
-
3
γ
2
(
1
-
ν
)
+
2
(
2
+
ν
2
)
-
3
γ
(
1
+
ν
)
-
3
γ
2
(
1
-
ν
)
-
2
(
2
+
ν
2
)
-
3
b
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
-
2
(
4
-
ν
2
)
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
+
ν
)
-
3
(
1
-
ν
)
+
2
γ
2
(
2
+
ν
2
)
-
3
γ
(
1
-
3
ν
)
-
3
(
1
-
ν
)
-
2
γ
2
(
2
+
ν
2
)
3
γ
(
1
+
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)
3
γ
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
-
2
(
4
-
ν
2
)
3
γ
(
1
+
ν
)
-
3
γ
2
(
1
-
ν
)
-
2
(
2
+
ν
2
)
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
3
γ
(
1
+
ν
)
3
(
1
-
ν
)
-
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
-
3
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)
3
γ
(
1
-
3
ν
)
-
3
γ
2
(
1
-
ν
)
+
2
(
2
+
ν
2
)
3
(
1
-
ν
)
+
2
γ
2
(
4
-
ν
2
)
-
3
γ
(
1
+
ν
)
3
γ
2
(
1
-
ν
)
+
2
(
4
-
ν
2
)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 whoseith mass M
i and stiffness K
i matrices are(24)
M
i
=
ρ
1
H
1
+
ρ
2
H
2
+
ρ
3
H
3
b
a
i
420
·
156
22
a
i
54
-
13
a
i
22
a
i
4
a
i
2
13
a
i
-
3
a
i
2
54
13
a
i
156
-
22
a
i
-
13
a
i
-
3
a
i
2
-
22
a
i
4
a
i
2
,
(25)
K
i
=
B
a
i
3
12
6
a
i
-
12
6
a
i
6
a
i
4
a
i
2
-
6
a
i
2
a
i
2
-
12
-
6
a
i
12
-
6
a
i
6
a
i
2
a
i
2
-
6
a
i
4
a
i
2
,respectively, where a
i is the length of the ith finite element. The complex flexural stiffness B of (25) is given by (2) B
=
B
eq
∗ for the RKU model and by (19) B
=
B
K
∗ for the new thick beam model. The discretisation is also made with 60 finite elements along span.


## 3.2. Extraction of Eigenvalues

The equation from which the complex eigenvalues of the system under study can be obtained is(26)
-
λ
r
∗
M
+
K
∗
ω
r
ϕ
r
∗
=
0
,where λ
r
∗ and ϕ
r
∗ are the complex eigenvalue and eigenvector of the rth mode, respectively, M is the mass matrix, and K
∗ is the complex stiffness matrix, which is dependent on frequency. This ω
r is the real part of the square root of the complex eigenvalue λ
r
∗:(27)
ω
r
=
Re
λ
r
∗
,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(28)
-
λ
r
,
0
M
+
K
∗
0
ϕ
r
,
0
=
0
,and then the undamped eigensolutions λ
r
,
0 and ϕ
r
,
0 can be obtained. Then, Nelson’s method [35] is used in order to calculate the eigenvector derivative ϕ
r
′. From this eigenvector derivative and taking into account the variation of the complex stiffness for the obtained eigenfrequency,(29)
Δ
K
∗
ω
r
,
0
=
K
∗
ω
r
,
0
-
K
∗
(
0
)
,and by means of Taylor’s series approach, a complex finite increment Δ
ϕ
r
∗ of the eigenvector is obtained. The complex eigenvector can be approximated with(30)
ϕ
r
∗
=
ϕ
r
,
0
+
Δ
ϕ
r
∗
,with which the complex eigenvalue λ
r
∗ is estimated according to(31)
λ
r
∗
=
ϕ
r
∗
H
Κ
∗
(
ω
r
)
ϕ
r
∗
ϕ
r
∗
H
M
ϕ
r
∗
,where (
·
)
H denotes the Hermitian transpose operator, that is, the complex conjugate transposition. Equations (29)–(31) can be iterated making use of the new eigenfrequency ω
r 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ω
-
r derived from (27) and loss factor derived from(32)
λ
r
∗
=
ω
-
r
2
1
+
i
η
r
.can be computed. The corresponding results for the three thicknesses and for the three models under study are shown in Tables 2–4.Table 2
Modal properties of the sandwich cantilever beam withH
2
=
1 mm.
ω
-
1 (rad/s)
η
1
ω
-
2 (rad/s)
η
2
ω
-
3 (rad/s)
η
3
2D model
411.1
0.024
2523
0.049
6851
0.047
Present model
412.2
0.019
2547
0.037
6912
0.043
RKU model
412.4
0.018
2552
0.034
6952
0.037Table 3
Modal properties of the sandwich cantilever beam withH
2
=
5 mm.
ω
-
1 (rad/s)
η
1
ω
-
2 (rad/s)
η
2
ω
-
3 (rad/s)
η
3
2D model
1003
0.080
5659
0.113
13987
0.091
Present model
1013
0.066
5909
0.094
14504
0.103
RKU model
1019
0.052
6051
0.068
15461
0.063Table 4
Modal properties of the sandwich cantilever beam withH
2
=
10 mm.
ω
-
1 (rad/s)
η
1
ω
-
2 (rad/s)
η
2
ω
-
3 (rad/s)
η
3
2D model
1561
0.118
8106
0.138
18828
0.099
Present model
1589
0.098
8714
0.127
19779
0.131
RKU model
1613
0.069
9291
0.076
22875
0.064It can be pointed out that the results for the natural frequencyω
r 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 Table3), 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η
r, 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 η
r depends on frequency according to (21), and the natural frequencies for the models are not the same. Instead of modal loss factor η
r, the amplitudes of the resonance peaks will be compared in the next section.


## 3.3. Analysis of the Frequency Response Function

The matrix system for a cantilever beam is given by(33)
-
ω
k
2
M
+
K
∗
ω
k
U
k
∗
=
F
∗
,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
k
∗ at each desired frequency ω
k. 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 F
=
1
×
exp
⁡
(
i
ω
t
)
N is applied on the right-hand side of our system.Figure 4
Frequency response function for the free-end for the cantilever beam: (a)H
2
=
1
mm, (b) H
2
=
5
mm, and (c) H
2
=
10
mm.
(a)
(b)
(c)Figure4(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, Tables5–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.Table 5
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thicknessH
2
=
1 mm.
A
1
A
2
A
3
2D model
-
6.357
-
8.608
-
10.58
Present model
-
6.299
-
8.304
-
10.46
RKU model
-
6.203
-
8.211
-
10.31Table 6
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thicknessH
2
=
5 mm.
A
1
A
2
A
3
2D model
-
7.552
-
11.44
-
13.10
Present model
-
7.285
-
11.16
-
13.04
RKU model
-
7.058
-
10.87
-
12.68Table 7
Logarithm of the resonance peak amplitude (m) of the free edge with core layer thicknessH
2
=
10 mm.
A
1
A
2
A
3
2D model
-
9.120
-
12.71
-
14.15
Present model
-
8.812
-
12.48
-
14.13
RKU model
-
8.451
-
12.06
-
13.72Since 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.


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



---

*Source: 736256-2015-04-07.xml*