# Transverse Vibration of the Thin Plates: Frequency-Domain Spectral Element Modeling and Analysis

Authors: Ilwook Park; Usik Lee; Donghyun Park
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

It has been well known that exact closed-form solutions are not available for non-Levy-type plates. Thus, more accurate and efficient computational methods have been required for the plates subjected to arbitrary boundary conditions. This paper presents a frequency-domain spectral element model for the rectangular finite plate element. The spectral element model is developed by using two methods in combination: (1) the boundary splitting and (2) the super spectral element method in which the Kantorovich method-based finite strip element method and the frequency-domain waveguide method are utilized. The present spectral element model has nodes on four edges of the finite plate element, but no nodes inside. This can reduce the total number of degrees of freedom a lot to improve the computational efficiency significantly, when compared with the standard finite element method (FEM). The high solution accuracy and computational efficiency of the present spectral element model are evaluated by the comparison with exact solutions and the solutions by the standard FEM.


---

## Body

## 1. Introduction

Although the subject of dynamic analysis of the plate has a long history which spans more than two centuries, most theoretical solutions to this classical problem are considered to be at best approximate because it is difficult to obtain exact closed-form solutions which simultaneously satisfy the governing differential equations and the associated boundary conditions, except for the Levy-type plates [1, 2].The finite element method (FEM) is one of the most powerful computational methods to predict the dynamic responses or waves in a complex structure with arbitrary boundary conditions. It has been well known that the accuracy of FEM depends on the number of finite elements (or mesh sizes) employed in the analysis as the displacement fields in a finite element are represented by using simple polynomials which are independent of vibrating frequencies. As a high frequency excitation force can excite the higher vibration modes, the mesh sizes must be at least smaller than the wavelength of the highest vibration mode of interest in order to acquire sufficiently accurate dynamic responses at the high frequency. In such circumstances, the FEM will be computationally too expensive. Thus, the spectral element method (SEM), as a combination of the flexibility of the FEM and the accuracy of the continuum element methods, can be considered as an alternative solution technique.In the literature, there are two different element methods which have been confusingly called the same name “SEM.” The first one is the fast Fourier transform- (FFT-) based frequency-domain SEM [3, 4], where accurately formulated dynamic stiffness matrix is used as the finite element stiffness matrix. The second one is the time-domain SEM proposed by Patera [5], which is formulated in the time domain by using the Legendre or Chebyshev orthogonal polynomials as the shape functions in conjunction with the use of the Gauss-Lobatto-Legendre integration rule. Thus, the word “spectral” is time-wise for the frequency-domain SEM, while it is space-wise for the time-domain SEM. The spectral element approach considered in this study is the frequency-domain SEM and it will be called “SEM” throughout this paper.In the SEM, the exact dynamic stiffness matrix (often called spectral element matrix) formulated from free wave solutions that satisfy the governing differential equations in the frequency domain is used as the finite element stiffness matrix. Thus, the SEM can provide exact solutions by representing a whole uniform structure member as a single finite element, regardless of its size, and it has been well known as an exact element method. The assembly of finite elements to a global structure system can be conducted in an exactly analogous way as that used in the conventional FEM. Despite the outstanding features of the SEM, however, its applications have been limited mostly to one-dimensional (1D) structures [3, 4], the isotropic, orthotropic, or laminated composite Levy-type plates [6–12], the semi-infinite plates (e.g., bounded in the x-direction but unbounded in the y-direction) [13–17], the infinite plates [18], the plates with completely free edge boundary conditions [19], and the plates with prespecified boundary conditions on two opposite edges [20, 21]. The dynamic stiffness matrices or spectral element models introduced in the previous studies [6–21] are valid only for very limited plate elements with specific geometries or boundary conditions. Thus, they cannot be used as a general finite plate element to deal with any finite plates subjected to arbitrary boundary conditions.In the literature, there are very few spectral element models for two-dimensional (2D) structures such as the plates and membranes which have finite dimensions in bothx- and y-directions and are subjected to arbitrary boundary conditions. This is because, for such finite 2D structures, it is not easy to obtain exact free wave solutions in analytical forms which are required to formulate exact dynamic stiffness matrices or spectral element models. Thus, a few researchers [22–27] have introduced some approximate dynamic stiffness matrices or spectral element models for the transverse vibrations of plates. Kulla [22] introduced a dynamic stiffness matrix for a finite plate element, however without presenting the details of the final form of formulation in his paper. He presented the frequency-domain dynamic response in terms of a finite series of base functions which satisfy the governing differential equation, with approximately representing the continuous boundary conditions in terms of truncated Fourier series. Gupta and Meek [23] introduced the so-called dynamic element method, where the dynamic stiffness matrix for a finite plate element was formulated by expressing the shape functions in the truncated ascending power of natural frequency. Recently, Lee and Jang [24] introduced a spectral element model for a finite plate element by representing the displacement field in a finite plate element as a series of plane waves which are propagating in various directions.In this study, a new spectral element model is proposed for a finite rectangular plate element. The proposed spectral element model is formulated by using two methods in combination. The first one is the method of splitting [25] where the boundary conditions are split into several parts to form partial problems. The second method is the spectral super element method (SSEM) [20] in which the Kantorovich method-based finite strip element method (in one direction) and the SEM (in the other direction) are used in combination. The SSEM is used to obtain the solutions to each partial problem and the partial solutions are summed up to obtain the full solution to the original problem. The performance of the proposed spectral element model is evaluated by the comparison with exact solutions and the solutions by the standard FEM.


## 2. Spectral Element Model for a Finite Plate Element

The small amplitude transverse vibration of a thin plate can be represented by the equation of motion as follows:(1)
D
∂
4
w
∂
x
4
+
2
∂
4
w
∂
x
2
∂
y
2
+
∂
4
w
∂
y
4
+
ρ
∂
2
w
∂
t
2
=
f
x
,
y
,
t
,where w
(
x
,
y
,
t
) is the transverse displacement, f
(
x
,
y
,
t
) is the external force applied perpendicularly to the surface of the plate, ρ is the mass per unit area of the plate, and D
=
E
h
3
/
[
12
(
1
-
ν
2
)
] is the flexural bending rigidity of the plate, where E is the modulus of elasticity, ν is Poisson’s ratio, and h is the thickness of the plate. The boundary conditions associated with (1) are given by(2)
M
y
x
,
-
L
y
2
,
t
=
-
M
y
1
x
,
t
or
θ
y
x
,
-
L
y
2
,
t
=
θ
y
1
x
,
t
V
y
x
,
-
L
y
2
,
t
=
-
V
y
1
x
,
t
or
w
x
,
-
L
y
2
,
t
=
w
1
x
,
t
M
y
x
,
L
y
2
,
t
=
M
y
2
x
,
t
or
θ
y
x
,
L
y
2
,
t
=
θ
y
2
x
,
t
V
y
x
,
L
y
2
,
t
=
V
y
2
x
,
t
or
w
x
,
L
y
2
,
t
=
w
2
x
,
t
M
x
-
L
x
2
,
y
,
t
=
-
M
x
1
y
,
t
or
θ
x
-
L
x
2
,
y
,
t
=
θ
x
1
y
,
t
V
x
-
L
x
2
,
y
,
t
=
-
V
x
1
y
,
t
or
w
-
L
x
2
,
y
,
t
=
w
3
y
,
t
M
x
L
x
2
,
y
,
t
=
M
x
2
y
,
t
or
θ
x
L
x
2
,
y
,
t
=
θ
x
2
y
,
t
V
x
L
x
2
,
y
,
t
=
V
x
2
y
,
t
or
w
L
2
x
,
y
,
t
=
w
4
y
,
t
,where M
x and M
y are the resultant moments and V
x and V
y are the resultant transverse shear forces defined by(3)
M
x
=
-
D
∂
2
w
∂
x
2
+
ν
∂
2
w
∂
y
2
,
M
y
=
-
D
∂
2
w
∂
y
2
+
ν
∂
2
w
∂
x
2
,
V
x
=
-
D
∂
∂
x
∂
2
w
∂
x
2
+
ν
∂
2
w
∂
y
2
-
2
D
1
-
ν
∂
3
w
∂
x
∂
y
2
,
V
y
=
-
D
∂
∂
y
∂
2
w
∂
y
2
+
ν
∂
2
w
∂
x
2
-
2
D
1
-
ν
∂
3
w
∂
x
2
∂
yand θ
x and θ
y are the slopes defined by(4)
θ
x
=
∂
w
∂
x
,
θ
y
=
∂
w
∂
y
.To formulate spectral element model for a thin plate by following the general procedure introduced in [4], firstly all the time-domain quantities in the governing equations (1) and (2) are transformed into the frequency-domain quantities by using the discrete Fourier transform (DFT) theory. For instance, the transverse displacement w
(
x
,
y
,
t
) and the external force f
(
x
,
y
,
t
) can be represented in the spectral forms as(5)
w
x
,
y
,
t
,
f
x
,
y
,
t
=
1
M
∑
m
=
0
M
-
1
w
-
m
x
,
y
,
f
-
m
x
,
y
e
i
ω
m
t
,where i
=
-
1 represents the imaginary unit. In (5), M is the number of samples for the fast Fourier transforms (FFT) analysis and ω
m
=
2
π
m
/
T are the discrete frequencies up to the Nyquist frequency. In this paper, the overbarred quantities such as w
-
m
(
x
,
y
) and f
-
m
(
x
,
y
) represent spectral (Fourier) components of the corresponding time-domain quantities. In the following derivation, the overbars and the subscripts m will be omitted for the sake of brevity. By using the spectral form representation of all time-domain quantities, the time-domain equation of motion (1) can be transformed into the frequency domain as (6)
D
∂
4
w
∂
x
4
+
2
∂
4
w
∂
x
2
∂
y
2
+
∂
4
w
∂
y
4
-
ρ
ω
2
w
=
f
x
,
y
.Similarly, boundary conditions (2) can be also transformed in the frequency domain as(7)
M
y
x
,
-
L
y
2
=
-
M
y
1
x
or
θ
y
x
,
-
L
y
2
=
θ
y
1
x
V
y
x
,
-
L
y
2
=
-
V
y
1
x
or
w
x
,
-
L
y
2
=
w
1
x
M
y
x
,
L
y
2
=
M
y
2
x
or
θ
y
x
,
L
y
2
=
θ
y
2
x
V
y
x
,
L
y
2
=
V
y
2
x
or
w
x
,
L
y
2
=
w
2
x
M
x
-
L
x
2
,
y
=
-
M
x
1
y
or
θ
x
-
L
x
2
,
y
=
θ
x
1
y
V
x
-
L
x
2
,
y
=
-
V
x
1
y
or
w
-
L
x
2
,
y
=
w
3
y
M
x
L
x
2
,
y
=
M
x
2
y
or
θ
x
L
x
2
,
y
=
θ
x
2
y
V
x
L
x
2
,
y
=
V
x
2
y
or
w
L
2
x
,
y
=
w
4
y
.Secondly, we need to obtain free wave solution in the frequency domain for a rectangular thin plate element that has the dimensionsL
x and L
y in the x- and y-directions, respectively. To that end, we consider the homogeneous equation of motion that can be readily obtained, by removing the external force term f
(
x
,
y
) from (6), as follows:(8)
D
∂
4
w
∂
x
4
+
2
∂
4
w
∂
x
2
∂
y
2
+
∂
4
w
∂
y
4
-
ρ
ω
2
w
=
0
.The solution w
(
x
,
y
) to (8) can be approximately obtained by using the boundary splitting method [25] in conjunction with the use of SSEM [20].The boundary splitting method used in this study is illustrated in Figure1. The original problem shown in Figure 1(a) is represented by two partial problems, Problem A and Problem B, by splitting original boundary conditions as shown in Figures 1(b) and 1(c), respectively. In Figure 1, the geometric boundary conditions are presented in simplified forms by using the symbols defined by(9)
w
x
x
,
y
=
w
x
,
y
,
θ
x
x
,
y
=
∂
w
∂
x
T
,
w
y
x
,
y
=
w
x
,
y
,
θ
y
x
,
y
=
∂
w
∂
y
T
,
w
x
A
x
,
y
=
w
A
x
,
y
,
θ
x
A
x
,
y
=
∂
w
A
∂
x
T
,
w
y
A
x
,
y
=
w
A
x
,
y
,
θ
y
A
x
,
y
=
∂
w
A
∂
y
T
,
w
x
B
x
,
y
=
w
B
x
,
y
,
θ
x
B
x
,
y
=
∂
w
B
∂
x
T
,
w
y
B
x
,
y
=
w
B
x
,
y
,
θ
y
B
x
,
y
=
∂
w
B
∂
y
T
.Figure 1
Boundary splitting method used to derivew
(
x
,
y
) for a rectangular plate element subjected to arbitrary boundary conditions.ProblemA has the fixed boundary conditions on two parallel edges at y
=
-
L
y
/
2 and L
y
/
2 and its solution is represented by w
A
(
x
,
y
). On the other hand, Problem B has the fixed boundary conditions on two parallel edges at x
=
-
L
x
/
2 and L
x
/
2 and its solution is represented by w
B
(
x
,
y
). The solution w
(
x
,
y
) to the original problem can be then obtained by summing the solutions to Problem A and Problem B as follows:(10)
w
x
,
y
=
w
A
x
,
y
+
w
B
x
,
y
.For the derivations ofw
A
(
x
,
y
) and w
B
(
x
,
y
) in the next subsections, the weak form of (8) can be obtained as(11)
∫
x
∫
y
D
∂
2
w
∂
x
2
+
ν
∂
2
w
∂
y
2
δ
∂
2
w
∂
x
2
+
2
1
-
ν
D
∂
2
w
∂
x
∂
y
δ
∂
2
w
∂
x
∂
y
+
D
∂
2
w
∂
y
2
+
ν
∂
2
w
∂
x
2
δ
∂
2
w
∂
y
2
-
ρ
ω
2
w
δ
w
d
x
d
y
=
0
.


### 2.1. Derivation ofw
A
(
x
,
y
)

To obtain the solutionw
A
(
x
,
y
) for Problem A by using the SSEM, the rectangular finite plate element that has the dimensions L
x and L
y in the x- and y-directions, respectively, is divided into a total of N
y finite strip elements in the y-direction as shown in Figure 2(a). The eth finite strip element that has the width of l
y
(
e
)
=
y
e
-
y
e
-
1 in the y-direction is shown in Figure 2(b).Figure 2
Finite strip element representation of a rectangular finite plate element subjected to null boundary conditions aty
=
-
L
y
/
2 and L
y
/
2, where ● represents nodes.
(a)
(b)The displacement field in theeth finite strip element can be represented by(12)
w
A
e
x
,
y
=
Y
A
e
y
u
A
e
x
y
e
-
1
≤
y
≤
y
e
,where Y
A
(
e
)
(
y
) is the one-by-four interpolation function matrix and u
A
(
e
)
(
x
) are the nodal line degrees of freedom (DOFs) defined by(13)
Y
A
e
y
=
Y
A
1
e
y
,
Y
A
2
e
y
,
u
A
e
x
=
w
y
A
e
-
1
x
w
y
A
e
x
,where(14)
Y
A
1
e
y
=
l
y
e
-
3
y
-
y
e
2
2
y
+
y
e
-
3
y
e
-
1
,
l
y
e
-
2
y
-
y
e
2
y
-
y
e
-
1
,
Y
A
2
e
y
=
-
l
y
e
-
3
y
-
y
e
-
1
2
2
y
+
y
e
-
1
-
3
y
e
,
l
y
e
-
2
y
-
y
e
-
1
2
y
-
y
e
,
w
y
A
e
x
=
w
y
A
x
,
y
e
=
w
A
e
x
,
y
e
,
θ
y
A
e
x
,
y
e
T
,
w
y
A
e
-
1
x
=
w
y
A
x
,
y
e
-
1
=
w
A
e
x
,
y
e
-
1
,
θ
y
A
e
x
,
y
e
-
1
T
.By using (12), the displacement field w
A
(
x
,
y
) in the whole domain of the finite plate element can be represented as(15)
w
A
x
,
y
=
Y
A
y
u
A
x
-
L
y
2
≤
y
≤
L
y
2
,where(16)
Y
A
y
=
L
A
1
y
,
L
A
2
y
,
…
,
L
A
e
y
,
…
,
L
A
N
y
-
2
y
,
L
A
N
y
-
1
y
,
u
A
x
=
w
y
A
1
x
,
w
y
A
2
x
,
…
,
w
y
A
e
x
,
…
,
w
y
A
N
y
-
2
x
,
w
y
A
N
y
-
1
x
T
.In (16), the following definition is used:(17)
L
A
e
y
=
h
A
e
y
Y
A
2
e
y
+
h
A
e
+
1
y
Y
A
1
e
+
1
y
e
=
1,2
,
3
,
…
,
N
y
-
1
,where (18)
h
A
e
y
=
H
y
-
y
e
-
1
-
H
y
-
y
e
,where H
(
y
) is the Heaviside unit step function. For (15), the null boundary conditions at y
=
-
L
y
/
2 and y
=
L
y
/
2 have been applied.By substituting (15) into (11), we can obtain(19)
A
A
4
∂
4
u
A
x
∂
x
4
+
A
A
2
∂
2
u
A
x
∂
x
2
+
A
A
0
-
ω
2
M
A
u
A
x
=
0
,where(20)
A
A
0
=
D
Λ
y
4
,
A
A
2
=
ν
D
Λ
y
3
+
Λ
y
3
T
-
2
1
-
ν
D
Λ
y
2
,
A
A
4
=
D
Λ
y
1
,
M
A
=
ρ
Λ
y
1
,
Λ
y
1
=
∫
-
L
y
/
2
L
y
/
2
Y
A
T
Y
A
d
y
,
Λ
y
2
=
∫
-
L
y
/
2
L
y
/
2
∂
Y
A
T
∂
y
∂
Y
A
∂
y
d
y
,
Λ
y
3
=
∫
-
L
y
/
2
L
y
/
2
∂
Y
A
T
∂
y
Y
A
d
y
,
Λ
y
4
=
∫
-
L
y
/
2
L
y
/
2
∂
2
Y
A
T
∂
y
2
∂
2
Y
A
∂
y
2
d
y
.The constant matrices  Λ
y
1,  Λ
y
2,  Λ
y
3,   and  Λ
y
4  are provided in the Appendix.Assume that the solution to (19) can be written in the following form:(21)
u
A
x
=
1
r
2
⋮
r
2
N
y
-
1
a
e
k
x
x
-
1
/
2
k
-
x
L
x
=
r
a
e
k
x
x
-
1
/
2
k
-
x
L
x
,where a is constant, k
x is the wavenumber in the x-direction, and(22)
k
-
x
=
+
k
x
if
Re
⁡
k
x
>
0
-
k
x
if
Re
⁡
k
x
<
0
0
if
Re
⁡
k
x
=
0
.Substituting (21) into (19) gives an eigenvalue problem as follows:(23)
A
A
4
k
x
4
+
A
A
2
k
x
2
+
A
A
0
-
ω
2
M
A
r
A
=
0
.The dispersion relation (frequency-wavenumber relationship) can be obtained from (24)
det
⁡
A
A
4
k
x
4
+
A
A
2
k
x
2
+
A
A
0
-
ω
2
M
A
=
0
.By using the wavenumbersk
x
(
j
) (j
=
1,2
,
3
,
…
,
8
×
(
N
y
-
1
)) computed from (24), we can write the solution to (19) in the form as(25)
u
A
x
=
R
A
E
A
x
a
A
,where(26)
a
A
=
a
1
,
a
2
,
a
3
,
…
,
a
i
,
…
,
a
8
×
N
y
-
1
T
,
R
A
=
r
A
1
,
r
A
2
,
r
A
3
,
…
,
r
A
j
,
…
,
r
A
8
×
N
y
-
1
,
E
A
x
;
ω
=
diagonal
e
k
x
j
x
-
1
/
2
k
-
x
j
L
x
j
=
1,2
,
3
,
…
,
8
×
N
y
-
1
.In (26), r
A
(
j
) is the jth wave mode vector (or eigenvector) that can be readily computed from (23) by putting k
x
=
k
x
(
j
).The nodal DOFs at the nodes defined on the edges atx
=
-
L
x
/
2 and L
x
/
2 can be written in a vector form as(27)
d
A
=
d
A
1
d
A
2
=
u
A
-
L
x
2
∂
u
A
-
L
x
/
2
/
∂
x
_
u
A
L
x
/
2
∂
u
A
L
x
/
2
∂
x
,where(28)
u
A
-
L
x
2
=
w
A
1
1
,
w
A
1
2
,
…
,
w
A
1
e
,
…
,
w
A
1
N
y
-
2
,
w
A
1
N
y
-
1
T
,
∂
u
A
-
L
x
/
2
∂
x
=
θ
x
A
1
1
,
θ
x
A
1
2
,
…
,
θ
x
A
1
e
,
…
,
θ
x
A
1
N
y
-
2
,
θ
x
A
1
N
y
-
1
T
,
u
A
L
x
2
=
w
A
2
1
,
w
A
2
2
,
…
,
w
A
2
e
,
…
,
w
A
2
N
y
-
2
,
w
A
2
N
y
-
1
T
,
∂
u
A
L
x
/
2
∂
x
=
θ
x
A
2
1
,
θ
x
A
2
2
,
…
,
θ
x
A
2
e
,
…
,
θ
x
A
2
N
y
-
2
,
θ
x
A
2
N
y
-
1
T
,with(29)
w
A
1
e
=
w
y
A
e
-
L
x
2
=
w
A
e
-
L
x
2
,
y
e
,
θ
y
A
e
-
L
x
2
,
y
e
T
,
θ
x
A
1
e
=
∂
w
y
A
e
-
L
x
/
2
∂
x
=
θ
x
A
e
-
L
x
2
,
y
e
,
θ
x
y
A
e
-
L
x
2
,
y
e
T
,
w
A
2
e
=
w
y
A
e
L
x
2
=
w
A
e
L
x
2
,
y
e
,
∂
w
A
e
L
x
/
2
,
y
e
∂
y
T
,
θ
x
A
2
e
=
∂
w
y
A
e
L
x
/
2
∂
x
=
θ
x
A
e
L
x
2
,
y
e
,
θ
x
y
A
e
L
x
2
,
y
e
T
.By substituting (25) into (27), the nodal DOFs vector d
A can be written in terms of the constant vector a
A as follows:(30)
d
A
=
R
A
E
A
-
L
x
2
R
A
K
A
E
A
-
L
x
2
R
A
E
A
+
L
x
2
R
A
K
A
E
A
+
L
x
2
a
A
=
H
A
ω
a
A
.The constant vector a
A can be removed from (25) by using (29) to obtain(31)
u
A
x
=
X
A
x
;
ω
d
A
,where(32)
X
A
x
;
ω
=
R
A
E
A
x
H
A
ω
-
1
.Finally, substituting (31) into (15), the solution w
A
(
x
,
y
) can be obtained in the form as(33)
w
A
x
,
y
=
N
A
x
,
y
;
ω
d
A
,where(34)
N
A
x
,
y
;
ω
=
Y
A
y
X
A
x
;
ω
.


### 2.2. Derivation ofw
B
(
x
,
y
)

As ProblemB can be obtained from Problem A by simply rotating the coordinates (
x
,
y
) 90° clockwise, we can find the solution w
B
(
x
,
y
) for Problem B by using the same procedure as that used in the previous subsection for Problem A. The only differences from Problem A are as follows: (1) the fixed boundary conditions are placed at x
=
-
L
x
/
2 and L
x
/
2, and (2) the finite plate element is divided into N
x finite strip elements in the x-direction (see Figure 3). Thus, without repeating similar procedure of solution derivation, the solution w
B
(
x
,
y
) for Problem B can be written as follows:(35)
w
B
x
,
y
=
X
B
x
Y
B
y
;
ω
d
B
=
N
B
x
,
y
;
ω
d
B
,where(36)
d
B
=
d
B
1
∣
d
B
2
T
=
w
B
1
1
,
w
B
1
2
,
…
,
w
B
1
e
,
…
,
w
B
1
N
x
-
2
,
w
B
1
N
x
-
1
,
θ
y
B
1
1
,
θ
y
B
1
2
,
…
,
θ
y
B
1
e
,
…
,
⁡
θ
y
B
1
N
x
-
2
,
θ
y
B
1
N
x
-
1
∣
w
B
2
1
,
w
B
2
2
,
…
,
w
B
2
e
,
…
,
w
B
2
N
y
-
2
,
w
B
2
N
y
-
1
,
θ
y
B
2
1
,
θ
y
B
2
2
,
…
,
θ
y
B
2
e
,
…
,
θ
y
B
2
N
x
-
2
,
θ
y
B
2
N
x
-
1
T
,
X
B
x
=
L
B
1
x
,
L
B
2
x
,
…
,
L
B
e
x
,
…
,
L
B
N
x
-
2
x
,
L
B
N
x
-
1
x
,
Y
B
y
;
ω
=
R
B
E
B
y
H
B
ω
-
1
.Figure 3
Finite strip element representation of a rectangular finite plate element subjected to null boundary conditions aty
=
-
L
x
/
2 and L
x
/
2, where ● represents nodes.
(a)
(b)In (36), the following definitions are used:(37)
w
B
1
e
=
w
x
B
e
-
L
y
2
=
w
B
e
x
e
,
-
L
y
2
,
θ
x
B
e
x
e
,
-
L
y
2
T
,
θ
y
B
1
e
=
∂
w
x
B
e
-
L
x
/
2
∂
x
=
θ
y
B
e
x
e
,
-
L
y
2
,
θ
x
y
B
e
x
e
,
-
L
y
2
T
,
w
B
2
e
=
w
x
B
e
L
x
2
=
w
B
e
x
e
,
L
y
2
,
θ
x
B
e
x
e
,
L
y
2
T
,
θ
y
B
2
e
=
∂
w
x
B
e
L
x
/
2
∂
x
=
θ
y
B
e
x
e
,
L
y
2
,
θ
x
y
B
e
x
e
,
L
y
2
T
,
(38)
L
B
e
y
=
h
B
e
y
X
B
2
e
x
+
h
B
e
+
1
y
X
B
1
e
+
1
x
,
R
B
=
r
B
1
,
r
B
2
,
r
B
3
,
…
,
r
B
j
,
…
,
r
B
8
N
x
-
1
,
E
B
y
;
ω
=
diagonal
e
k
y
j
y
-
1
/
2
k
-
y
j
L
y
j
=
1,2
,
3
,
…
,
8
×
N
x
-
1
,
H
B
ω
=
R
B
E
B
-
L
y
2
R
B
K
B
E
B
-
L
y
2
R
B
E
B
L
y
2
R
B
K
B
E
B
L
y
2
,where(39)
X
B
1
e
=
l
x
e
-
3
x
-
x
e
2
2
x
-
3
x
e
-
1
+
x
e
,
l
x
e
-
2
x
-
x
e
-
1
x
-
x
e
2
,
X
B
2
e
=
-
l
x
e
-
3
x
-
x
e
-
1
2
2
x
+
x
e
-
1
-
3
x
e
,
l
x
e
-
2
x
-
x
e
x
-
x
e
-
1
2
,
h
B
e
x
=
H
x
-
x
e
-
1
-
H
x
-
x
e
,
k
-
y
j
=
+
k
y
j
if
Re
⁡
k
y
j
>
0
-
k
y
j
if
Re
⁡
k
y
j
<
0
0
if
Re
⁡
k
y
j
=
0
.Thejth wavenumber k
y
(
j
) and the corresponding wave mode vector r
B
(
j
), which appear in (38), can be computed by solving the eigenvalue problem given by (40)
A
B
4
k
y
j
4
+
A
B
2
k
y
j
2
+
A
B
0
-
ω
2
M
B
r
B
j
=
0
,where(41)
A
B
0
=
D
Λ
x
4
,
A
B
2
=
ν
D
Λ
x
3
+
Λ
x
3
T
-
2
1
-
ν
D
Λ
x
2
,
A
B
4
=
D
Λ
x
1
,
M
B
=
ρ
Λ
x
1
.The constant matrices  Λ
x
1,  Λ
x
2,  Λ
x
3,   and  Λ
x
4  are provided in the Appendix.


### 2.3. Derivation ofw
(
x
,
y
)

The solutionw
(
x
,
y
) for the original problem shown in Figure 1(a) can be obtained by substituting (33) and (35) into (10) as follows:(42)
w
x
,
y
=
N
A
x
,
y
;
ω
d
A
+
N
B
x
,
y
;
ω
d
B
=
N
x
,
y
;
ω
d
,where  d
(
ω
) is the 8
×
(
N
x
+
N
y
-
2
)-by-one spectral nodal DOFs vector defined by (43)
d
=
d
A
d
Band N
(
x
,
y
;
ω
) is the one-by-8
×
(
N
x
+
N
y
-
2
) dynamic shape function matrix defined by(44)
N
x
,
y
;
ω
=
N
A
x
,
y
;
ω
N
B
x
,
y
;
ω
.


### 2.4. Formulation of Spectral Element Equation

The weak form of the governing equation of forced vibration, (6), can be derived as follows:(45)
∫
x
∫
y
D
∂
2
w
∂
x
2
+
ν
∂
2
w
∂
y
2
δ
∂
2
w
∂
x
2
+
2
1
-
ν
D
∂
2
w
∂
x
∂
y
δ
∂
2
w
∂
x
∂
y
+
D
∂
2
w
∂
y
2
+
ν
∂
2
w
∂
x
2
δ
∂
2
w
∂
y
2
-
ρ
ω
2
w
δ
w
d
x
d
y
=
∫
x
∫
y
f
x
,
y
δ
w
d
x
d
y
+
∫
y
V
x
1
δ
w
-
L
x
2
,
y
d
y
+
∫
y
V
x
2
δ
w
L
x
2
,
y
d
y
+
∫
x
V
y
1
δ
w
x
,
-
L
y
2
d
x
+
∫
x
V
y
2
δ
w
x
,
L
y
2
d
x
+
∫
y
M
x
1
δ
∂
w
-
L
x
/
2
,
y
∂
x
d
y
+
∫
y
M
x
2
δ
∂
w
L
x
/
2
,
y
∂
x
d
y
+
∫
x
M
y
1
δ
∂
w
x
,
-
L
y
/
2
∂
y
d
x
+
∫
x
M
y
2
δ
∂
w
x
,
L
y
/
2
∂
y
d
x
,where V
x
1
(
y
), V
x
2
(
y
), V
y
1
(
x
), and V
y
2
(
x
) are the transverse shear forces applied on four boundary edges. Similarly, M
x
1
(
y
), M
x
2
(
y
), M
y
1
(
x
), and M
y
2
(
x
) are the bending moments applied on four boundary edges.By substituting (42) into (44), we can obtain the spectral element equation in the following form:(46)
S
ω
d
ω
=
f
1
ω
+
f
2
ω
,where (47)
S
ω
=
H
-
T
ω
D
1
ω
+
D
2
ω
+
D
3
ω
+
D
4
ω
+
D
4
T
ω
-
ω
2
D
ρ
ω
H
-
1
ω
=
H
-
T
ω
D
ω
H
-
1
ω
,
(48)
f
1
ω
=
∫
x
∫
y
f
x
,
y
N
T
x
,
y
d
x
d
y
,
f
2
ω
=
∫
y
V
x
1
N
T
-
1
2
L
x
,
y
d
y
+
∫
y
V
x
2
N
T
1
2
L
x
,
y
d
y
+
∫
x
V
y
1
N
T
x
,
-
1
2
L
y
d
x
+
∫
x
V
y
2
N
T
x
,
1
2
L
y
d
x
+
∫
y
M
x
1
∂
N
T
-
1
/
2
L
x
,
y
∂
x
d
y
+
∫
y
M
x
2
∂
N
T
1
/
2
L
x
,
y
∂
x
d
y
+
∫
x
M
y
1
∂
N
T
x
,
-
1
/
2
L
y
∂
y
d
x
+
∫
x
M
y
2
∂
N
T
x
,
1
/
2
L
y
∂
y
d
x
.
8
(
N
x
+
N
y
-
2
)-by-8
(
N
x
+
N
y
-
2
) matrix S
(
ω
) is the frequency dependent dynamic stiffness matrix, which is often called spectral element matrix. The symbols used in (47) are defined by(49)
H
ω
=
H
A
ω
0
0
H
B
ω
,
(50)
D
1
=
D
K
A
2
R
A
T
Λ
A
1
R
A
K
A
2
.
∗
Γ
1
Γ
4
R
B
.
∗
K
A
2
R
A
T
Γ
6
T
symmetric
R
B
T
Λ
A
4
R
B
.
∗
Γ
5
,
D
2
=
D
R
A
T
Λ
A
4
R
A
.
∗
Γ
1
Γ
2
R
B
K
B
2
.
∗
R
A
T
Γ
8
T
symmetric
K
B
2
R
B
T
Λ
B
1
R
B
K
B
2
.
∗
Γ
1
,
D
3
=
2
1
-
ν
D
K
A
R
A
T
Λ
A
2
R
A
K
A
.
∗
Γ
1
Γ
3
R
B
K
B
.
∗
K
A
R
A
T
Γ
7
T
symmetric
K
B
R
B
T
Λ
B
2
R
B
K
B
.
∗
Γ
5
,
D
4
=
ν
D
K
A
2
R
A
T
Λ
A
3
R
A
.
∗
Γ
1
Γ
2
R
B
K
B
2
.
∗
K
A
2
R
A
T
Γ
6
T
Γ
8
R
A
.
∗
R
B
T
Γ
4
T
R
B
T
Λ
B
3
T
R
B
K
B
2
.
∗
Γ
5
,
D
ρ
=
ρ
R
A
T
Λ
A
1
R
A
.
∗
Γ
1
Γ
2
R
B
.
∗
R
A
T
Γ
6
T
symmetric
R
B
T
Λ
B
1
R
B
.
∗
Γ
5
, where(51)
Γ
1
=
∫
x
vect
⁡
E
A
vect
⁡
E
A
T
d
x
,
Γ
2
=
∫
x
vect
⁡
E
A
X
B
d
x
,
Γ
3
=
∫
x
vect
⁡
E
A
∂
X
B
∂
x
d
x
,
Γ
4
=
∫
x
vect
⁡
E
A
∂
2
X
B
∂
x
2
d
x
,
Γ
5
=
∫
y
vect
⁡
E
B
vect
⁡
E
B
T
d
y
,
Γ
6
=
∫
y
vect
⁡
E
B
Y
A
d
y
,
Γ
7
=
∫
y
vect
⁡
E
B
∂
Y
A
∂
y
d
y
,
Γ
8
=
∫
y
vect
⁡
E
B
∂
2
Y
A
∂
y
2
d
y
.In (50), the symbol (.
∗) denotes the element-wise matrix multiplication defined in Matlab [28] as follows:(52)
A
.
∗
B
=
C
,where the components of matrix C are defined by(53)
C
i
j
=
A
i
j
×
B
i
j
.


## 2.1. Derivation ofw
A
(
x
,
y
)

To obtain the solutionw
A
(
x
,
y
) for Problem A by using the SSEM, the rectangular finite plate element that has the dimensions L
x and L
y in the x- and y-directions, respectively, is divided into a total of N
y finite strip elements in the y-direction as shown in Figure 2(a). The eth finite strip element that has the width of l
y
(
e
)
=
y
e
-
y
e
-
1 in the y-direction is shown in Figure 2(b).Figure 2
Finite strip element representation of a rectangular finite plate element subjected to null boundary conditions aty
=
-
L
y
/
2 and L
y
/
2, where ● represents nodes.
(a)
(b)The displacement field in theeth finite strip element can be represented by(12)
w
A
e
x
,
y
=
Y
A
e
y
u
A
e
x
y
e
-
1
≤
y
≤
y
e
,where Y
A
(
e
)
(
y
) is the one-by-four interpolation function matrix and u
A
(
e
)
(
x
) are the nodal line degrees of freedom (DOFs) defined by(13)
Y
A
e
y
=
Y
A
1
e
y
,
Y
A
2
e
y
,
u
A
e
x
=
w
y
A
e
-
1
x
w
y
A
e
x
,where(14)
Y
A
1
e
y
=
l
y
e
-
3
y
-
y
e
2
2
y
+
y
e
-
3
y
e
-
1
,
l
y
e
-
2
y
-
y
e
2
y
-
y
e
-
1
,
Y
A
2
e
y
=
-
l
y
e
-
3
y
-
y
e
-
1
2
2
y
+
y
e
-
1
-
3
y
e
,
l
y
e
-
2
y
-
y
e
-
1
2
y
-
y
e
,
w
y
A
e
x
=
w
y
A
x
,
y
e
=
w
A
e
x
,
y
e
,
θ
y
A
e
x
,
y
e
T
,
w
y
A
e
-
1
x
=
w
y
A
x
,
y
e
-
1
=
w
A
e
x
,
y
e
-
1
,
θ
y
A
e
x
,
y
e
-
1
T
.By using (12), the displacement field w
A
(
x
,
y
) in the whole domain of the finite plate element can be represented as(15)
w
A
x
,
y
=
Y
A
y
u
A
x
-
L
y
2
≤
y
≤
L
y
2
,where(16)
Y
A
y
=
L
A
1
y
,
L
A
2
y
,
…
,
L
A
e
y
,
…
,
L
A
N
y
-
2
y
,
L
A
N
y
-
1
y
,
u
A
x
=
w
y
A
1
x
,
w
y
A
2
x
,
…
,
w
y
A
e
x
,
…
,
w
y
A
N
y
-
2
x
,
w
y
A
N
y
-
1
x
T
.In (16), the following definition is used:(17)
L
A
e
y
=
h
A
e
y
Y
A
2
e
y
+
h
A
e
+
1
y
Y
A
1
e
+
1
y
e
=
1,2
,
3
,
…
,
N
y
-
1
,where (18)
h
A
e
y
=
H
y
-
y
e
-
1
-
H
y
-
y
e
,where H
(
y
) is the Heaviside unit step function. For (15), the null boundary conditions at y
=
-
L
y
/
2 and y
=
L
y
/
2 have been applied.By substituting (15) into (11), we can obtain(19)
A
A
4
∂
4
u
A
x
∂
x
4
+
A
A
2
∂
2
u
A
x
∂
x
2
+
A
A
0
-
ω
2
M
A
u
A
x
=
0
,where(20)
A
A
0
=
D
Λ
y
4
,
A
A
2
=
ν
D
Λ
y
3
+
Λ
y
3
T
-
2
1
-
ν
D
Λ
y
2
,
A
A
4
=
D
Λ
y
1
,
M
A
=
ρ
Λ
y
1
,
Λ
y
1
=
∫
-
L
y
/
2
L
y
/
2
Y
A
T
Y
A
d
y
,
Λ
y
2
=
∫
-
L
y
/
2
L
y
/
2
∂
Y
A
T
∂
y
∂
Y
A
∂
y
d
y
,
Λ
y
3
=
∫
-
L
y
/
2
L
y
/
2
∂
Y
A
T
∂
y
Y
A
d
y
,
Λ
y
4
=
∫
-
L
y
/
2
L
y
/
2
∂
2
Y
A
T
∂
y
2
∂
2
Y
A
∂
y
2
d
y
.The constant matrices  Λ
y
1,  Λ
y
2,  Λ
y
3,   and  Λ
y
4  are provided in the Appendix.Assume that the solution to (19) can be written in the following form:(21)
u
A
x
=
1
r
2
⋮
r
2
N
y
-
1
a
e
k
x
x
-
1
/
2
k
-
x
L
x
=
r
a
e
k
x
x
-
1
/
2
k
-
x
L
x
,where a is constant, k
x is the wavenumber in the x-direction, and(22)
k
-
x
=
+
k
x
if
Re
⁡
k
x
>
0
-
k
x
if
Re
⁡
k
x
<
0
0
if
Re
⁡
k
x
=
0
.Substituting (21) into (19) gives an eigenvalue problem as follows:(23)
A
A
4
k
x
4
+
A
A
2
k
x
2
+
A
A
0
-
ω
2
M
A
r
A
=
0
.The dispersion relation (frequency-wavenumber relationship) can be obtained from (24)
det
⁡
A
A
4
k
x
4
+
A
A
2
k
x
2
+
A
A
0
-
ω
2
M
A
=
0
.By using the wavenumbersk
x
(
j
) (j
=
1,2
,
3
,
…
,
8
×
(
N
y
-
1
)) computed from (24), we can write the solution to (19) in the form as(25)
u
A
x
=
R
A
E
A
x
a
A
,where(26)
a
A
=
a
1
,
a
2
,
a
3
,
…
,
a
i
,
…
,
a
8
×
N
y
-
1
T
,
R
A
=
r
A
1
,
r
A
2
,
r
A
3
,
…
,
r
A
j
,
…
,
r
A
8
×
N
y
-
1
,
E
A
x
;
ω
=
diagonal
e
k
x
j
x
-
1
/
2
k
-
x
j
L
x
j
=
1,2
,
3
,
…
,
8
×
N
y
-
1
.In (26), r
A
(
j
) is the jth wave mode vector (or eigenvector) that can be readily computed from (23) by putting k
x
=
k
x
(
j
).The nodal DOFs at the nodes defined on the edges atx
=
-
L
x
/
2 and L
x
/
2 can be written in a vector form as(27)
d
A
=
d
A
1
d
A
2
=
u
A
-
L
x
2
∂
u
A
-
L
x
/
2
/
∂
x
_
u
A
L
x
/
2
∂
u
A
L
x
/
2
∂
x
,where(28)
u
A
-
L
x
2
=
w
A
1
1
,
w
A
1
2
,
…
,
w
A
1
e
,
…
,
w
A
1
N
y
-
2
,
w
A
1
N
y
-
1
T
,
∂
u
A
-
L
x
/
2
∂
x
=
θ
x
A
1
1
,
θ
x
A
1
2
,
…
,
θ
x
A
1
e
,
…
,
θ
x
A
1
N
y
-
2
,
θ
x
A
1
N
y
-
1
T
,
u
A
L
x
2
=
w
A
2
1
,
w
A
2
2
,
…
,
w
A
2
e
,
…
,
w
A
2
N
y
-
2
,
w
A
2
N
y
-
1
T
,
∂
u
A
L
x
/
2
∂
x
=
θ
x
A
2
1
,
θ
x
A
2
2
,
…
,
θ
x
A
2
e
,
…
,
θ
x
A
2
N
y
-
2
,
θ
x
A
2
N
y
-
1
T
,with(29)
w
A
1
e
=
w
y
A
e
-
L
x
2
=
w
A
e
-
L
x
2
,
y
e
,
θ
y
A
e
-
L
x
2
,
y
e
T
,
θ
x
A
1
e
=
∂
w
y
A
e
-
L
x
/
2
∂
x
=
θ
x
A
e
-
L
x
2
,
y
e
,
θ
x
y
A
e
-
L
x
2
,
y
e
T
,
w
A
2
e
=
w
y
A
e
L
x
2
=
w
A
e
L
x
2
,
y
e
,
∂
w
A
e
L
x
/
2
,
y
e
∂
y
T
,
θ
x
A
2
e
=
∂
w
y
A
e
L
x
/
2
∂
x
=
θ
x
A
e
L
x
2
,
y
e
,
θ
x
y
A
e
L
x
2
,
y
e
T
.By substituting (25) into (27), the nodal DOFs vector d
A can be written in terms of the constant vector a
A as follows:(30)
d
A
=
R
A
E
A
-
L
x
2
R
A
K
A
E
A
-
L
x
2
R
A
E
A
+
L
x
2
R
A
K
A
E
A
+
L
x
2
a
A
=
H
A
ω
a
A
.The constant vector a
A can be removed from (25) by using (29) to obtain(31)
u
A
x
=
X
A
x
;
ω
d
A
,where(32)
X
A
x
;
ω
=
R
A
E
A
x
H
A
ω
-
1
.Finally, substituting (31) into (15), the solution w
A
(
x
,
y
) can be obtained in the form as(33)
w
A
x
,
y
=
N
A
x
,
y
;
ω
d
A
,where(34)
N
A
x
,
y
;
ω
=
Y
A
y
X
A
x
;
ω
.


## 2.2. Derivation ofw
B
(
x
,
y
)

As ProblemB can be obtained from Problem A by simply rotating the coordinates (
x
,
y
) 90° clockwise, we can find the solution w
B
(
x
,
y
) for Problem B by using the same procedure as that used in the previous subsection for Problem A. The only differences from Problem A are as follows: (1) the fixed boundary conditions are placed at x
=
-
L
x
/
2 and L
x
/
2, and (2) the finite plate element is divided into N
x finite strip elements in the x-direction (see Figure 3). Thus, without repeating similar procedure of solution derivation, the solution w
B
(
x
,
y
) for Problem B can be written as follows:(35)
w
B
x
,
y
=
X
B
x
Y
B
y
;
ω
d
B
=
N
B
x
,
y
;
ω
d
B
,where(36)
d
B
=
d
B
1
∣
d
B
2
T
=
w
B
1
1
,
w
B
1
2
,
…
,
w
B
1
e
,
…
,
w
B
1
N
x
-
2
,
w
B
1
N
x
-
1
,
θ
y
B
1
1
,
θ
y
B
1
2
,
…
,
θ
y
B
1
e
,
…
,
⁡
θ
y
B
1
N
x
-
2
,
θ
y
B
1
N
x
-
1
∣
w
B
2
1
,
w
B
2
2
,
…
,
w
B
2
e
,
…
,
w
B
2
N
y
-
2
,
w
B
2
N
y
-
1
,
θ
y
B
2
1
,
θ
y
B
2
2
,
…
,
θ
y
B
2
e
,
…
,
θ
y
B
2
N
x
-
2
,
θ
y
B
2
N
x
-
1
T
,
X
B
x
=
L
B
1
x
,
L
B
2
x
,
…
,
L
B
e
x
,
…
,
L
B
N
x
-
2
x
,
L
B
N
x
-
1
x
,
Y
B
y
;
ω
=
R
B
E
B
y
H
B
ω
-
1
.Figure 3
Finite strip element representation of a rectangular finite plate element subjected to null boundary conditions aty
=
-
L
x
/
2 and L
x
/
2, where ● represents nodes.
(a)
(b)In (36), the following definitions are used:(37)
w
B
1
e
=
w
x
B
e
-
L
y
2
=
w
B
e
x
e
,
-
L
y
2
,
θ
x
B
e
x
e
,
-
L
y
2
T
,
θ
y
B
1
e
=
∂
w
x
B
e
-
L
x
/
2
∂
x
=
θ
y
B
e
x
e
,
-
L
y
2
,
θ
x
y
B
e
x
e
,
-
L
y
2
T
,
w
B
2
e
=
w
x
B
e
L
x
2
=
w
B
e
x
e
,
L
y
2
,
θ
x
B
e
x
e
,
L
y
2
T
,
θ
y
B
2
e
=
∂
w
x
B
e
L
x
/
2
∂
x
=
θ
y
B
e
x
e
,
L
y
2
,
θ
x
y
B
e
x
e
,
L
y
2
T
,
(38)
L
B
e
y
=
h
B
e
y
X
B
2
e
x
+
h
B
e
+
1
y
X
B
1
e
+
1
x
,
R
B
=
r
B
1
,
r
B
2
,
r
B
3
,
…
,
r
B
j
,
…
,
r
B
8
N
x
-
1
,
E
B
y
;
ω
=
diagonal
e
k
y
j
y
-
1
/
2
k
-
y
j
L
y
j
=
1,2
,
3
,
…
,
8
×
N
x
-
1
,
H
B
ω
=
R
B
E
B
-
L
y
2
R
B
K
B
E
B
-
L
y
2
R
B
E
B
L
y
2
R
B
K
B
E
B
L
y
2
,where(39)
X
B
1
e
=
l
x
e
-
3
x
-
x
e
2
2
x
-
3
x
e
-
1
+
x
e
,
l
x
e
-
2
x
-
x
e
-
1
x
-
x
e
2
,
X
B
2
e
=
-
l
x
e
-
3
x
-
x
e
-
1
2
2
x
+
x
e
-
1
-
3
x
e
,
l
x
e
-
2
x
-
x
e
x
-
x
e
-
1
2
,
h
B
e
x
=
H
x
-
x
e
-
1
-
H
x
-
x
e
,
k
-
y
j
=
+
k
y
j
if
Re
⁡
k
y
j
>
0
-
k
y
j
if
Re
⁡
k
y
j
<
0
0
if
Re
⁡
k
y
j
=
0
.Thejth wavenumber k
y
(
j
) and the corresponding wave mode vector r
B
(
j
), which appear in (38), can be computed by solving the eigenvalue problem given by (40)
A
B
4
k
y
j
4
+
A
B
2
k
y
j
2
+
A
B
0
-
ω
2
M
B
r
B
j
=
0
,where(41)
A
B
0
=
D
Λ
x
4
,
A
B
2
=
ν
D
Λ
x
3
+
Λ
x
3
T
-
2
1
-
ν
D
Λ
x
2
,
A
B
4
=
D
Λ
x
1
,
M
B
=
ρ
Λ
x
1
.The constant matrices  Λ
x
1,  Λ
x
2,  Λ
x
3,   and  Λ
x
4  are provided in the Appendix.


## 2.3. Derivation ofw
(
x
,
y
)

The solutionw
(
x
,
y
) for the original problem shown in Figure 1(a) can be obtained by substituting (33) and (35) into (10) as follows:(42)
w
x
,
y
=
N
A
x
,
y
;
ω
d
A
+
N
B
x
,
y
;
ω
d
B
=
N
x
,
y
;
ω
d
,where  d
(
ω
) is the 8
×
(
N
x
+
N
y
-
2
)-by-one spectral nodal DOFs vector defined by (43)
d
=
d
A
d
Band N
(
x
,
y
;
ω
) is the one-by-8
×
(
N
x
+
N
y
-
2
) dynamic shape function matrix defined by(44)
N
x
,
y
;
ω
=
N
A
x
,
y
;
ω
N
B
x
,
y
;
ω
.


## 2.4. Formulation of Spectral Element Equation

The weak form of the governing equation of forced vibration, (6), can be derived as follows:(45)
∫
x
∫
y
D
∂
2
w
∂
x
2
+
ν
∂
2
w
∂
y
2
δ
∂
2
w
∂
x
2
+
2
1
-
ν
D
∂
2
w
∂
x
∂
y
δ
∂
2
w
∂
x
∂
y
+
D
∂
2
w
∂
y
2
+
ν
∂
2
w
∂
x
2
δ
∂
2
w
∂
y
2
-
ρ
ω
2
w
δ
w
d
x
d
y
=
∫
x
∫
y
f
x
,
y
δ
w
d
x
d
y
+
∫
y
V
x
1
δ
w
-
L
x
2
,
y
d
y
+
∫
y
V
x
2
δ
w
L
x
2
,
y
d
y
+
∫
x
V
y
1
δ
w
x
,
-
L
y
2
d
x
+
∫
x
V
y
2
δ
w
x
,
L
y
2
d
x
+
∫
y
M
x
1
δ
∂
w
-
L
x
/
2
,
y
∂
x
d
y
+
∫
y
M
x
2
δ
∂
w
L
x
/
2
,
y
∂
x
d
y
+
∫
x
M
y
1
δ
∂
w
x
,
-
L
y
/
2
∂
y
d
x
+
∫
x
M
y
2
δ
∂
w
x
,
L
y
/
2
∂
y
d
x
,where V
x
1
(
y
), V
x
2
(
y
), V
y
1
(
x
), and V
y
2
(
x
) are the transverse shear forces applied on four boundary edges. Similarly, M
x
1
(
y
), M
x
2
(
y
), M
y
1
(
x
), and M
y
2
(
x
) are the bending moments applied on four boundary edges.By substituting (42) into (44), we can obtain the spectral element equation in the following form:(46)
S
ω
d
ω
=
f
1
ω
+
f
2
ω
,where (47)
S
ω
=
H
-
T
ω
D
1
ω
+
D
2
ω
+
D
3
ω
+
D
4
ω
+
D
4
T
ω
-
ω
2
D
ρ
ω
H
-
1
ω
=
H
-
T
ω
D
ω
H
-
1
ω
,
(48)
f
1
ω
=
∫
x
∫
y
f
x
,
y
N
T
x
,
y
d
x
d
y
,
f
2
ω
=
∫
y
V
x
1
N
T
-
1
2
L
x
,
y
d
y
+
∫
y
V
x
2
N
T
1
2
L
x
,
y
d
y
+
∫
x
V
y
1
N
T
x
,
-
1
2
L
y
d
x
+
∫
x
V
y
2
N
T
x
,
1
2
L
y
d
x
+
∫
y
M
x
1
∂
N
T
-
1
/
2
L
x
,
y
∂
x
d
y
+
∫
y
M
x
2
∂
N
T
1
/
2
L
x
,
y
∂
x
d
y
+
∫
x
M
y
1
∂
N
T
x
,
-
1
/
2
L
y
∂
y
d
x
+
∫
x
M
y
2
∂
N
T
x
,
1
/
2
L
y
∂
y
d
x
.
8
(
N
x
+
N
y
-
2
)-by-8
(
N
x
+
N
y
-
2
) matrix S
(
ω
) is the frequency dependent dynamic stiffness matrix, which is often called spectral element matrix. The symbols used in (47) are defined by(49)
H
ω
=
H
A
ω
0
0
H
B
ω
,
(50)
D
1
=
D
K
A
2
R
A
T
Λ
A
1
R
A
K
A
2
.
∗
Γ
1
Γ
4
R
B
.
∗
K
A
2
R
A
T
Γ
6
T
symmetric
R
B
T
Λ
A
4
R
B
.
∗
Γ
5
,
D
2
=
D
R
A
T
Λ
A
4
R
A
.
∗
Γ
1
Γ
2
R
B
K
B
2
.
∗
R
A
T
Γ
8
T
symmetric
K
B
2
R
B
T
Λ
B
1
R
B
K
B
2
.
∗
Γ
1
,
D
3
=
2
1
-
ν
D
K
A
R
A
T
Λ
A
2
R
A
K
A
.
∗
Γ
1
Γ
3
R
B
K
B
.
∗
K
A
R
A
T
Γ
7
T
symmetric
K
B
R
B
T
Λ
B
2
R
B
K
B
.
∗
Γ
5
,
D
4
=
ν
D
K
A
2
R
A
T
Λ
A
3
R
A
.
∗
Γ
1
Γ
2
R
B
K
B
2
.
∗
K
A
2
R
A
T
Γ
6
T
Γ
8
R
A
.
∗
R
B
T
Γ
4
T
R
B
T
Λ
B
3
T
R
B
K
B
2
.
∗
Γ
5
,
D
ρ
=
ρ
R
A
T
Λ
A
1
R
A
.
∗
Γ
1
Γ
2
R
B
.
∗
R
A
T
Γ
6
T
symmetric
R
B
T
Λ
B
1
R
B
.
∗
Γ
5
, where(51)
Γ
1
=
∫
x
vect
⁡
E
A
vect
⁡
E
A
T
d
x
,
Γ
2
=
∫
x
vect
⁡
E
A
X
B
d
x
,
Γ
3
=
∫
x
vect
⁡
E
A
∂
X
B
∂
x
d
x
,
Γ
4
=
∫
x
vect
⁡
E
A
∂
2
X
B
∂
x
2
d
x
,
Γ
5
=
∫
y
vect
⁡
E
B
vect
⁡
E
B
T
d
y
,
Γ
6
=
∫
y
vect
⁡
E
B
Y
A
d
y
,
Γ
7
=
∫
y
vect
⁡
E
B
∂
Y
A
∂
y
d
y
,
Γ
8
=
∫
y
vect
⁡
E
B
∂
2
Y
A
∂
y
2
d
y
.In (50), the symbol (.
∗) denotes the element-wise matrix multiplication defined in Matlab [28] as follows:(52)
A
.
∗
B
=
C
,where the components of matrix C are defined by(53)
C
i
j
=
A
i
j
×
B
i
j
.


## 3. Numerical Results and Discussion

To evaluate the performance of the present spectral element model for the finite thin plate elements, two types of thin plates are considered as the numerical examples: (1) Example 1 (a simply supported uniform square plate as shown in Figure4(a)) and (2) Example 2 (a clamped square plate with a square cutout inside, all edges of which are also clamped as shown in Figure 5(a)). We assumed that two example plates are made of the same material (aluminum) and their material properties are given as follows: Young’s modulus E
=
69 GPa, Poisson’s ratio ν
=
0.33, and mass density ρ
=
2700 kg/m3. The thickness of the plates is given by h
=
0.001 m and the other dimensions are shown in Figures 4 and 5.Figure 4
Example 1: (a) simply supported square plate, (b) FEM model, and (c) SEM models, where ● represents nodes (units: m).
(a)
Example 1
(b)
FEM model
(c)
SEM modelFigure 5
Example 2: (a) clamped square plate with a square cutout inside, (b) FEM model, and (c) SEM model, where ● represents nodes (units: m).
(a)
Example 2
(b)
FEM model
(c)
SEM modelEvaluation of the present spectral element method (denoted by “SEM”) has been conducted by comparing the natural frequencies obtained by the SEM with those obtained by the standard finite element method (denoted by “FEM”). In this study, the four-node 12-DOF conforming rectangular finite element model [27] has been used for the FEM results. Exact solutions are available in analytical forms for the Levy-type plates as shown in Figure 4. The exact natural frequencies of a simply supported square plate are given as follows [26]:(54)
f
m
n
=
π
m
2
+
n
2
2
L
2
D
ρ
Hz
m
,
n
=
1,2
,
3
,
…
,where L is the dimension of the square plate and (
m
,
n
) indicates the mode number. Thus, exact natural frequencies given by (54) (denoted by “Exact”) were also used to evaluate the accuracies of the present SEM and the standard FEM.Table1 compares the natural frequencies of the simply supported square plate (i.e., Example 1) obtained by the present SEM with those by the standard FEM as well as with exact natural frequencies given by (54). For the SEM results, the square plate shown in Figure 4(a) was represented by the one-element model as shown in Figure 4(c) and the number of finite strip elements in the x
- and y-directions (i.e., N
x and N
y) was increased until the natural frequency of mode (
4,4
) was sufficiently converged to exact solution. Similarly, for the FEM results, the number of finite elements used in the standard finite element analysis was also increased until the natural frequency of mode (
4,4
) was sufficiently converged to exact solution. Table 1 shows that the FEM results certainly converge to exact solutions when the total number of DOFs used in the analysis is increased up to about 43000. Meanwhile, the SEM results are found to converge to exact solutions when the total number of DOFs employed in the one-element model is increased up to about 400. This investigation shows that the same level of solution accuracy achieved by the standard FEM can be achieved by the present SEM by using an extremely small number of DOFs.Table 1
Comparison of the natural frequencies (Hz) of the simply supported square plate (Example1) obtained by the exact theory, FEM [27], and SEM.
Modes
Exact
FEM
SEM
n
E = 50 × 50
n
D = 7399
100 × 10029799
120 × 12042959
n
E = 1
n
D = 72
1152
1392
(1, 1)
4.857
4.857
4.857
4.866
4.866
4.857
4.857
(1, 2)
12.14
12.14
12.14
12.16
12.16
12.14
12.14
(2, 2)
19.43
19.42
19.43
19.46
19.46
19.43
19.43
(1, 3)
24.28
24.26
24.28
24.28
24.28
4.857
24.28
(2, 3)
31.57
31.54
31.57
31.60
31.60
31.58
31.57
(4, 1)
41.28
41.20
41.26
41.30
41.30
41.29
41.28
(3, 3)
43.71
43.68
43.71
43.75
43.75
43.73
43.71
(4, 2)
48.57
48.49
48.52
48.53
48.53
48.60
48.57
(4, 3)
60.71
60.59
60.64
60.72
60.72
60.73
60.71
(4, 4)
77.71
77.51
77.60
77.66
77.66
77.73
77.71
Note.
n
E: total number of finite elements; n
D: total number of degrees of freedom.Exact solutions are not available in analytical forms for Example 2 (i.e., the clamped square plate with a square cutout inside as shown in Figure5(a)). Thus, for the evaluation of the present SEM, the natural frequencies obtained by the present SEM are compared in Table 2 with the results obtained by the standard FEM and a commercial finite element analysis package ANSYS [29]. For the SEM results, the example plate shown in Figure 5(a) was represented by eight same-sized square finite elements as shown in Figure 5(c) and the number of finite strip elements in the x- and y-directions (i.e., N
x and N
y) for each square finite element was increased until the eighth natural frequency was sufficiently converged to exact solution. For both ANSYS and FEM results, the numbers of finite elements used in the analysis were also increased until the eighth natural frequency is sufficiently converged to exact solution. The natural frequencies by ANSYS and the standard FEM are found to converge very close to each other when the total number of DOFs employed for ANSYS and FEM analyses was increased up to about 78800 and 39400, respectively. On the other hand, Table 2 shows that the same level of solution accuracy can be achieved by the present SEM by using an extremely small number of DOFs (i.e., about 1900). Figure 6 shows the lowest six mode shapes obtained by the present SEM.Table 2
Comparison of the natural frequencies (Hz) of the clamped square plate with a square cutout inside (Example2) obtained by ANSYS, FEM [27], and SEM.
Modes
ANSYS
FEM
SEM
n
E = 12801
n
D = 78726
n
E = 3201
n
D = 9603
720122323
1280139363
n
E = 8
n
D = 608
8928
81,888
1
55.13
55.10
55.13
55.13
55.15
55.13
55.13
2
55.73
55.69
55.73
55.73
55.75
55.73
55.73
3
56.62
56.55
56.57
56.62
56.67
56.64
56.62
4
65.52
65.48
65.49
65.51
65.59
65.54
65.53
5
69.60
69.53
69.57
69.60
69.68
69.63
69.60
6
77.19
77.05
77.15
77.18
77.30
77.24
77.20
7
83.18
83.05
83.15
83.17
83.33
83.28
83.20
8
94.63
94.43
94.56
94.60
94.71
94.65
94.64
Note.
n
E: total number of finite elements; n
D: total number of degrees of freedom.Figure 6
The lowest six mode shapes of the clamped square plate with a square cutout inside (Example 2).Lastly, we considered a clamped uniform square plate of dimensionL
=
1 m (without a square cutout inside) and assumed that point impulse loading was applied at the center of the square plate. Figure 7 shows the impulse-induced waves predicted by the present SEM at times t
=
500 μs, 1000 μs, and 1500 μs. The SEM results shown in Figure 7 are obtained by using a total of 297 DOFs. Figure 7 shows that the impulse-induced waves are generated at the middle part of the square plate and they propagate radially in a circular wave form until they reach the nearest boundary point of the plate. After reaching the nearest boundary point, the wave forms become complex more and more due to the superposition of the forward propagating waves and the backward propagating waves that are continuously reflected at all boundary edges.Figure 7
The impulse-induced waves in the clamped uniform square plate predicted at times (a)t
=
500 μs, (b) t
=
1000 μs, and (c) t
=
1500 μs.
(a)
(b)
(c)Figure8 compares the impulse-induced transient dynamic responses predicted by the present SEM at three locations, x
/
(
L
/
2
)
=
0.2, 0.22, and 0.24, on the x-axis of the square plate. Figure 8 shows that the wave form measured at x
/
(
L
/
2
)
=
0.2 propagates toward next two locations, x
/
(
L
/
2
)
=
0.22 and x
/
(
L
/
2
)
=
0.24, with time.Figure 8
Comparison of the impulse-induced transient dynamic responses predicted at three locations on thex-axis of the clamped uniform square plate.Figure9 compares the impulse-induced transient dynamic response predicted by the present SEM with exact solution for a uniform infinite plate given by [30](55)
w
r
,
t
=
1
2
π
ρ
ρ
D
π
2
-
Si
r
2
4
t
ρ
D
,where Si
(
x
)
=
∫
0
x
(
sin
⁡
z
)
z
-
1
d
z is the sine integral and r is the distance from the point at which a unit impulse is applied.Figure 9
Comparison of the impulse-induced transient dynamic responses (at a locationx
/
(
L
/
2
)
=
0.2 on the x-axis of the clamped uniform square plate) predicted by the present SEM with exact analytical solution for an infinite plate [30].The impulse-induced transient dynamic responses were predicted at a locationx
/
(
L
/
2
)
=
0.2 on the x-axis of the uniform plate. As the impulse-induced waves propagate radially in a circular wave form until they reach the nearest boundary point of the plate, the transient dynamic response predicted by the present SEM at the early time of excitation is found to be almost identical to that predicted by exact analytical solution for a uniform infinite plate [30]. However, Figure 9 shows that the transient dynamic response predicted by the present SEM deviates significantly from the exact analytical solution for a uniform infinite plate [30] from the time about t
=
720 μs when the waves reflected at the plate edges propagate backward to be superposed at the response measurement location x
/
(
L
/
2
)
=
0.2.Figure10 shows the comparison of the vibration responses obtained by the present SEM and the standard FEM [27]. The vibration responses were predicted at two locations x
/
(
L
/
2
)
=
0.2 and 0.4 on the x-axis of the clamped uniform square plate. The SEM and FEM results were obtained by using a total of 117 and 1083 DOFs, respectively. The CPU times required for the SEM-based and FEM-based dynamic simulations were 258 seconds and 2602 seconds, respectively. For both dynamic simulations, a standard desktop PC that is equipped with two sockets of Intel Xeon E5-2620 v2 processors and 256 GB of DDR3 RAM memory clocked at 1600 MHz was used. From this investigation, we can conclude that the SEM-based dynamic simulation is in general much faster than the FEM-based dynamic simulation, for instance, about ten times faster for the example problem considered in this study.Figure 10
Comparison of the impulse-induced vibration responses predicted by present SEM and FEM [27]: (a) at x
/
(
L
/
2
)
=
0.2 and (b) at x
/
(
L
/
2
)
=
0.4 on the x-axis of the clamped uniform square plate.
(a)
(b)


## 4. Conclusions

In this paper, a frequency-domain spectral element model for finite thin plates is proposed. The proposed spectral element model is developed by using two methods in combination: the boundary splitting and the super spectral element method (SSEM). In SSEM, the Kantorovich method-based finite strip element method and the frequency-domain waveguide method (or the one-dimensional spectral element modeling) are utilized. As a result, the proposed spectral element model has nodes on four edges of the finite plate element, but no nodes inside. This can result in a significant reduction of the total number of degrees of freedom used in the proposed SEM, when compared with the standard FEM, to result in the significant improvement of the computational efficiency. Some example problems are considered to demonstrate extremely high solution accuracy and computational efficiency of the proposed SEM by the comparison with the solutions by various solution techniques: exact solutions available in the literature, the standard FEM, and a commercial finite element analysis package ANSYS.



---

*Source: 541276-2015-12-14.xml*