# A Hybrid Finite Element-Fourier Spectral Method for Vibration Analysis of Structures with Elastic Boundary Conditions

Authors: Wan-You Li; Binglin Lv; Huajiang Ouyang; Jingtao Du; Haijun Zhou; Donghua Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A novel hybrid method, which simultaneously possesses the efficiency of Fourier spectral method (FSM) and the applicability of the finite element method (FEM), is presented for the vibration analysis of structures with elastic boundary conditions. The FSM, as one type of analytical approaches with excellent convergence and accuracy, is mainly limited to problems with relatively regular geometry. The purpose of  the current study is to extend the FSM to problems with irregular geometry via the FEM and attempt to take full advantage of the FSM and the conventional FEM for structural vibration problems. The computational domain of general shape is divided into several subdomains firstly, some of which are represented by the FSM while the rest by the FEM. Then, fictitious springs are introduced for connecting these subdomains. Sufficient details are given to describe the development of such a hybrid method. Numerical examples of a one-dimensional Euler-Bernoulli beam and a two-dimensional rectangular plate show that the present method has good accuracy and efficiency. Further, one irregular-shaped plate which consists of one rectangular plate and one semi-circular plate also demonstrates the capability of the present method applied to irregular structures.


---

## Body

## 1. Introduction

The needs for engineers to accurately predict performance of dynamic systems and to produce optimal designs as well as fast progress in hardware performance and constant decrease in the price of computers all lead to increasing popularity and application of the FEM in engineering. Though the FEM is applicable to problems of any geometric and material properties and boundary conditions, it can be very expensive and time consuming to solve large scale finite element models using a desktop computer, especially when the frequency involved is relatively high. When a parametric study is required in the analysis of some problems, such as for the rotor shaft-oil film bearing-flexible foundation system including several different components, remeshing and reanalysis are needed each time [1–3]. Additionally, the FEM model would become impractically large if wave propagation characteristics of a structure at higher frequencies are sought [4].Besides simplifying a model to reduce computational workload, many attempts have been made within the last decades in order to replace or improve the conventional FEM with other methods [5–19]. Some of these are able to solve mid-frequency vibration problems with some degrees of success. Apparently, techniques for solving mid-frequency vibration problems of general boundary conditions at high accuracy are still needed [20]. How to obtain an accurate model with fewer degrees of freedom is one of the main issues in structural dynamic analysis. Recently, Fourier spectral method (FSM) which was put forward by Li [21] has been widely applied to various vibration analyses of regular structures such as straight beams [21], multispan beams [22], and in-plane [23] and out-of-plane [24, 25] vibration of rectangular plates. At the early stage of this method, the generalized Fourier series are first substituted into the governing differential equations with relevant boundary equations which are derived from the force equilibrium relationship, and then the Galerkin discretization procedure is incorporated to transform the overall equations into a standard matrix eigenvalue problem. And recently, one more general formulation based on an energy principle where all the unknown coefficients are regarded as generalized independent coordinates [26] has been introduced, using FSM, to vibration analysis of, for example, single rectangular plates with arbitrary nonuniform elastic edge restraints [27], and more complex structures such as elastically coupled rectangular plates [28] and circular cylindrical shells [29]. FSM is even applied to vibroacoustic analysis of a three-dimensional acoustic cavity with both the displacement of the plate and the sound pressure constructed in the form of Fourier series with supplemental terms [30]. As a novel Fourier expansion method, the convergence speed of the FSM is improved substantially compared with the traditional Fourier series method. It is a suitable systematic method for solving a wide range of engineering and scientific problems. However, compared with local methods such as the FEM, as an analytical method which uses global basis functions, the FSM also shows intrinsic difficulties in satisfying boundary conditions for irregular-shaped domains like other global methods [16] such as spectral methods or spectral finite element methods [13, 14].For many problems such as ships and aeroplanes, the domain considered is geometrically complex and the application of the available spectral techniques to such structures would become difficult or even impossible. The main purpose of this paper is to combine the advantages of the FSM and FEM to develop a new hybrid method for problems with partly regular-shaped boundaries. Firstly, the structure is divided into two zones: one with regular boundaries described by the FSM and the other with irregular boundaries discretised by the FEM. At their interface, the coupling between these two zones is implemented by introducing several fictitious springs on the common edges. These two subsystems are coupled together using Hamilton’s equation. After these two methods are combined, it can be found that the computational efficiency is enhanced compared with the standard FEM, and limitations that the FSM can hardly be applied to complex structures are also overcome. Sections2 and 3 give, respectively, the theoretical formulations and numerical examples of one-dimensional (beam) and two-dimensional (plate) structures based on the new method.


## 2. Theoretical Formulations of 1-D Structures

### 2.1. Governing Equation of Beam Vibration

Consider an Euler-Bernoulli beam elastically supported at its both ends. The equations of motion as well as the corresponding boundary conditions can be derived based on the energy principle, in which the total potential energyV and total kinetic energy T are written as
(1)
V
=
1
2
∫
0
L
FSM
D
(
d
2
w
d
x
2
)
2
d
x
+
(
(
d
w
d
x
)
2
1
2
k
w
0
w
2
+
1
2
K
w
0
(
d
w
d
x
)
2
)
|
x
=
0
+
(
1
2
k
w
1
w
2
+
1
2
K
w
1
(
d
w
d
x
)
2
)
|
x
=
L
FSM
,
T
=
1
2
ω
2
ρ
A
∫
0
L
FSM
w
2
d
x
,
where L
FSM is the beam length, w is the flexural displacement, ω is the angular frequency, and D, ρ, and A are, respectively, the flexural rigidity, the mass density, and the cross-sectional area; k
0 and k
1 are the linear spring constants, and K
0 and K
1 are the rotational spring constants at x
=
0 and x
=
L
FSM, respectively. If FEM is used, the equation for free vibration can be written as
(2)
(
[
K
FEM
]
-
ω
2
[
M
FEM
]
)
Λ
=
{
0
}
,
in which Λ
=
{
w
1
,
θ
1
,
w
2
,
θ
2
,
…
,
w
S
,
θ
S
}
T and w
i and θ
i are, respectively, the deflection and rotation at the ith node.For a beam with elastic boundary conditions, it is not easy to determine its eigenfunctions which can be used as displacement functions. In this study, the displacement function for FSM will be sought in the form of a simple series expansion as(3)
w
(
x
)
=
∑
m
=
0
∞
a
m
cos
⁡
λ
L
m
x
+
∑
l
=
1
4
(
c
l
ξ
l
(
x
)
)
,
where λ
L
m
=
m
π
/
L
FSM. When cosine functions are used, the subsequent mathematical derivation becomes much easier. However, the elastic boundary conditions cannot be satisfied by such a series and discontinuities in the first and third derivatives occur at the boundaries. Four supplemental functions ξ
l are introduced to remove the potential discontinuities at both ends; namely,
(4)
ξ
1
(
x
)
=
9
L
FSM
4
π
sin
(
π
x
2
L
FSM
)
-
L
FSM
12
π
sin
(
3
π
x
2
L
FSM
)
,
ξ
2
(
x
)
=
-
9
L
FSM
4
π
cos
⁡
(
π
x
2
L
FSM
)
-
L
FSM
12
π
cos
⁡
(
3
π
x
2
L
FSM
)
,
ξ
3
(
x
)
=
L
FSM
3
π
3
sin
(
π
x
2
L
FSM
)
-
L
FSM
3
3
π
3
sin
(
3
π
x
2
L
FSM
)
,
ξ
4
(
x
)
=
-
L
FSM
3
π
3
cos
⁡
(
π
x
2
L
FSM
)
-
L
FSM
3
3
π
3
cos
⁡
(
3
π
x
2
L
FSM
)
.Actually, the choice of these auxiliary functions is not unique; however, the appropriate construction of these functions can simplify the subsequent mathematical formulation significantly. These functions can be chosen as long as the following conditions are satisfied:ξ
1
′
(
0
)
=
ξ
2
′
(
L
FSM
)
=
ξ
3
′′′
(
0
)
=
ξ
4
′′′
(
L
FSM
)
=
1 and all the rest first and third derivatives of these functions at these two end points are equal to zero. Detailed derivations and explanations can be found in [21, 27]. In order to find all the unknown coefficients in (3), which actually determines the structural dynamic behaviour, the Hamilton equation is employed here; namely,
(5)
δ
∫
t
0
t
1
(
T
-
V
)
d
t
=
0
,
where V denotes the total potential energy and T the total kinetic energy. Substituting (3) into the Hamilton equation and minimizing it with respect to each unknown Fourier coefficient and c
i will lead to the following equation in matrix form:
(6)
(
K
FSM
-
ω
2
M
FSM
)
Ψ
=
{
0
}
,
where Ψ
=
{
a
0
,
a
1
,
…
,
a
M
,
c
1
,
c
2
,
c
3
,
c
4
}
T and K
FSM and M
FSM are the system stiffness and mass matrices which are explicitly defined in Appendix A in supplementary material available online at http://dx.doi.org/10.1155/2014/834834. In (6), it is assumed that the series expansion (3) is truncated to m
=
M terms in practical calculation; namely, m
=
0,1
,
2
,
…
,
M.


### 2.2. Combination of the FEM and the FSM

After briefly reviewing the main principles of FEM and FSM, development of a hybrid method from these two will be described.As illustrated in Figure1, a beam structure is divided into two zones: the left is handled by the FEM and the right by the FSM, which are, respectively, named as FEM zone and FSM zone for convenience. To combine these two methods, two types of fictitious springs with linear and rotational spring coefficients of k
c and k
c are introduced at the interface where two zones meet. On the FEM side, it is assumed that the beam is divided into (
S
-
1
) elements, namely, with S nodes. The displacement and rotation at the interface are w
S
FEM and θ
S
FEM. On the other hand, the displacement and rotation at the interface from the FSM are
(7)
w
FSM
|
x
=
0
=
∑
m
=
0
∞
a
m
+
∑
l
=
1
4
(
c
l
ξ
l
(
0
)
)
,
θ
FSM
|
x
=
0
=
c
l
(
ξ
l
(
0
)
)
′
.Figure 1
Schematic of elastically restrained beam structure divided into FEM zone and FSM zone, with two types of fictitious springs introduced at zone interface.The potential energy stored in these two springs can now be calculated as(8)
V
w
c
=
1
2
k
c
(
w
S
FEM
-
w
FSM
|
x
=
0
)
2
,
(9)
V
θ
c
=
1
2
K
c
(
θ
S
FEM
-
θ
FSM
|
x
=
0
)
2
.When the spring constants of the fictitious springs involved are assigned very big values, the transverse displacements and rotations of the two subdomains at the interface will then be effectively equal.Substituting (8) and (9) into (5), the coupled matrix K
couple can then be obtained. The overall equation of motion can be derived as
(10)
(
[
K
FEM
0
0
K
FSM
]
+
[
K
couple
]
hhhh
-
ω
2
[
M
FEM
0
0
M
FSM
]
)
{
Λ
Ψ
}
=
{
0
}
,
where the elements of the coupled matrices are given in Appendix B.


### 2.3. Numerical Results and Discussion

In this section, several numerical examples are presented to show the accuracy of this hybrid method. Firstly, a straight beam with the following parameters is used in the calculation: beam length=
L
FEM
+
L
FSM
=
1
m , in which L
FEM and L
FSM represent the length of the FEM and FSM zones on the beam structure; Young’s modulus E
=
2.1
×
10
11
Pa; the second moment of inertia of the cross-section I
=
π
r
4
/
4 with r
=
0.01
m . The Fourier series are truncated to M
=
20 for the FSM, while the element number for the FEM is set as 50. In order to study the influence of the relative proportion of these two subdomains on the prediction results, the length of the FSM zone is varied from 0 to 1 with the increment of 0.2 while the overall length of the beam is kept unchanged as =
1
m . The results are given in Table 1 for a beam with clamped-clamped boundary condition, which can be readily achieved by setting all the end restraint stiffness into a very large number (the translational springs are k
w
0
=
5
×
1
0
7, k
w
1
=
5
×
1
0
7 while the rotation springs at the boundary are K
w
0
=
5
×
1
0
7, K
w
1
=
5
×
1
0
7 and the coupled springs are k
c
=
5
×
1
0
7, K
c
=
5
×
1
0
7) for FSM zone side and fixing the end node on FEM zone side. As can be observed, the relative proportion between these two computational zones almost has no effect on the accuracy of the predicted frequencies. Such good results indicate that the proposed hybrid method is effective.Table 1
Effect of the division ratio of the two computational domains on the modal frequencies of the beams structure with clamped-clamped boundary condition.
Mode order
μ
i
=
L
/
π
(
ω
i
ρ
A
/
EI
)
1
/
2, M = 20, E
n = 50
L
FSM (m)
Analytical results
1.0
0.8
0.6
0.4
0.2
0.0
L
FEM (m)
0.0
0.2
0.4
0.6
0.8
1.0
1
1.5056
1.5056
1.5056
1.5056
1.5056
1.5056
1.5055
2
2.4998
2.4998
2.4997
2.4998
2.4997
2.4998
2.4996
3
3.5000
3.5000
3.5000
3.5000
3.5000
3.5000
3.5000
4
4.5000
4.5000
4.5000
4.5000
4.5000
4.5000
4.4999
5
5.5000
5.5000
5.4999
5.4999
5.4999
5.4999
5.5000
6
6.5000
6.5000
6.4998
6.5000
6.5000
6.4998
6.5001
7
7.5000
7.5000
7.4999
7.4998
7.4999
7.4999
7.5002
8
8.5000
8.5000
8.5000
8.4996
8.4997
8.5001
8.5003
9
9.5000
9.5000
9.4999
9.5000
9.5000
9.5000
9.5004
10
10.5000
10.4999
10.4993
10.4992
10.4992
10.4995
10.5007Table2 compares the efficiency and accuracy of the FEM, the FSM, and the hybrid method. All these three methods are implemented on the same computational platform in MATLAB. A relative high frequency range up to the 20th mode is computed by these three approaches. Firstly, FEM and FSM are used to make prediction and it can be seen that FSM is much more efficient than FEM, since its dimension of final system matrix is much smaller than that of FEM. For the proposed hybrid method, it can be observed again that the prediction is accurate enough for capturing the modal characteristics of this free-free beam.Table 2
Comparison of computing time of the hybrid method with pure FEM and FEM with different truncated Fourier terms and/or FEM element numbers, in which a completely free-free beam structure is used.
Method
FSM truncatedterm, M
FEM elementnumber, E
n
The 20th modal frequency
Degrees of freedom
Computing time (s)
Analytical
—
—
19.50000000
—
—
FSM
80
—
19.50000044
85
0.016
FEM
—
600
19.50000073
1202
9.235
Hybrid
40
300
19.49999961
647
7.66
FSM
60
—
19.5000032
65
0.016
FEM
—
500
19.50000153
1002
5.133
Hybrid
30
250
19.50000005
537
3.947
FSM
50
—
19.5000114
55
0.015
FEM
—
250
19.50002436
502
0.64
Hybrid
25
125
19.50001615
282
0.499
FSM
30
—
19.50046727
35
<0.01
FEM
—
100
19.50094026
202
0.047
Hybrid
15
50
19.50069514
122
0.031In the first three lines of Table2, we see that in these three cases the accuracy class of the 20th modal frequency calculated is equal. The degrees of freedom of the matrix used for calculation are 85, 600, and 345 for FSM, FEM, and current hybrid method, respectively. Although the CPU times are not a good measure of an algorithm’s performance as they do not provide much insight about an algorithm’s behaviour and they are still presented here as references which are 0.016(s), 9.235(s), and 7.66(s), respectively, for the first three cases. Compared with FEM, since the FSM with rapid convergence and high efficiency is incorporated into this structure, the overall degrees of freedom will be reduced considerably, say, from 600 to 345 for the first three cases. The CPU time consumed is also reduced compared with pure FEM. The remaining cases in this table are all the same: the current method would reduce DOF significantly compared with pure FEM while the accuracy class of the frequency calculated is still kept in the same level.


## 2.1. Governing Equation of Beam Vibration

Consider an Euler-Bernoulli beam elastically supported at its both ends. The equations of motion as well as the corresponding boundary conditions can be derived based on the energy principle, in which the total potential energyV and total kinetic energy T are written as
(1)
V
=
1
2
∫
0
L
FSM
D
(
d
2
w
d
x
2
)
2
d
x
+
(
(
d
w
d
x
)
2
1
2
k
w
0
w
2
+
1
2
K
w
0
(
d
w
d
x
)
2
)
|
x
=
0
+
(
1
2
k
w
1
w
2
+
1
2
K
w
1
(
d
w
d
x
)
2
)
|
x
=
L
FSM
,
T
=
1
2
ω
2
ρ
A
∫
0
L
FSM
w
2
d
x
,
where L
FSM is the beam length, w is the flexural displacement, ω is the angular frequency, and D, ρ, and A are, respectively, the flexural rigidity, the mass density, and the cross-sectional area; k
0 and k
1 are the linear spring constants, and K
0 and K
1 are the rotational spring constants at x
=
0 and x
=
L
FSM, respectively. If FEM is used, the equation for free vibration can be written as
(2)
(
[
K
FEM
]
-
ω
2
[
M
FEM
]
)
Λ
=
{
0
}
,
in which Λ
=
{
w
1
,
θ
1
,
w
2
,
θ
2
,
…
,
w
S
,
θ
S
}
T and w
i and θ
i are, respectively, the deflection and rotation at the ith node.For a beam with elastic boundary conditions, it is not easy to determine its eigenfunctions which can be used as displacement functions. In this study, the displacement function for FSM will be sought in the form of a simple series expansion as(3)
w
(
x
)
=
∑
m
=
0
∞
a
m
cos
⁡
λ
L
m
x
+
∑
l
=
1
4
(
c
l
ξ
l
(
x
)
)
,
where λ
L
m
=
m
π
/
L
FSM. When cosine functions are used, the subsequent mathematical derivation becomes much easier. However, the elastic boundary conditions cannot be satisfied by such a series and discontinuities in the first and third derivatives occur at the boundaries. Four supplemental functions ξ
l are introduced to remove the potential discontinuities at both ends; namely,
(4)
ξ
1
(
x
)
=
9
L
FSM
4
π
sin
(
π
x
2
L
FSM
)
-
L
FSM
12
π
sin
(
3
π
x
2
L
FSM
)
,
ξ
2
(
x
)
=
-
9
L
FSM
4
π
cos
⁡
(
π
x
2
L
FSM
)
-
L
FSM
12
π
cos
⁡
(
3
π
x
2
L
FSM
)
,
ξ
3
(
x
)
=
L
FSM
3
π
3
sin
(
π
x
2
L
FSM
)
-
L
FSM
3
3
π
3
sin
(
3
π
x
2
L
FSM
)
,
ξ
4
(
x
)
=
-
L
FSM
3
π
3
cos
⁡
(
π
x
2
L
FSM
)
-
L
FSM
3
3
π
3
cos
⁡
(
3
π
x
2
L
FSM
)
.Actually, the choice of these auxiliary functions is not unique; however, the appropriate construction of these functions can simplify the subsequent mathematical formulation significantly. These functions can be chosen as long as the following conditions are satisfied:ξ
1
′
(
0
)
=
ξ
2
′
(
L
FSM
)
=
ξ
3
′′′
(
0
)
=
ξ
4
′′′
(
L
FSM
)
=
1 and all the rest first and third derivatives of these functions at these two end points are equal to zero. Detailed derivations and explanations can be found in [21, 27]. In order to find all the unknown coefficients in (3), which actually determines the structural dynamic behaviour, the Hamilton equation is employed here; namely,
(5)
δ
∫
t
0
t
1
(
T
-
V
)
d
t
=
0
,
where V denotes the total potential energy and T the total kinetic energy. Substituting (3) into the Hamilton equation and minimizing it with respect to each unknown Fourier coefficient and c
i will lead to the following equation in matrix form:
(6)
(
K
FSM
-
ω
2
M
FSM
)
Ψ
=
{
0
}
,
where Ψ
=
{
a
0
,
a
1
,
…
,
a
M
,
c
1
,
c
2
,
c
3
,
c
4
}
T and K
FSM and M
FSM are the system stiffness and mass matrices which are explicitly defined in Appendix A in supplementary material available online at http://dx.doi.org/10.1155/2014/834834. In (6), it is assumed that the series expansion (3) is truncated to m
=
M terms in practical calculation; namely, m
=
0,1
,
2
,
…
,
M.


## 2.2. Combination of the FEM and the FSM

After briefly reviewing the main principles of FEM and FSM, development of a hybrid method from these two will be described.As illustrated in Figure1, a beam structure is divided into two zones: the left is handled by the FEM and the right by the FSM, which are, respectively, named as FEM zone and FSM zone for convenience. To combine these two methods, two types of fictitious springs with linear and rotational spring coefficients of k
c and k
c are introduced at the interface where two zones meet. On the FEM side, it is assumed that the beam is divided into (
S
-
1
) elements, namely, with S nodes. The displacement and rotation at the interface are w
S
FEM and θ
S
FEM. On the other hand, the displacement and rotation at the interface from the FSM are
(7)
w
FSM
|
x
=
0
=
∑
m
=
0
∞
a
m
+
∑
l
=
1
4
(
c
l
ξ
l
(
0
)
)
,
θ
FSM
|
x
=
0
=
c
l
(
ξ
l
(
0
)
)
′
.Figure 1
Schematic of elastically restrained beam structure divided into FEM zone and FSM zone, with two types of fictitious springs introduced at zone interface.The potential energy stored in these two springs can now be calculated as(8)
V
w
c
=
1
2
k
c
(
w
S
FEM
-
w
FSM
|
x
=
0
)
2
,
(9)
V
θ
c
=
1
2
K
c
(
θ
S
FEM
-
θ
FSM
|
x
=
0
)
2
.When the spring constants of the fictitious springs involved are assigned very big values, the transverse displacements and rotations of the two subdomains at the interface will then be effectively equal.Substituting (8) and (9) into (5), the coupled matrix K
couple can then be obtained. The overall equation of motion can be derived as
(10)
(
[
K
FEM
0
0
K
FSM
]
+
[
K
couple
]
hhhh
-
ω
2
[
M
FEM
0
0
M
FSM
]
)
{
Λ
Ψ
}
=
{
0
}
,
where the elements of the coupled matrices are given in Appendix B.


## 2.3. Numerical Results and Discussion

In this section, several numerical examples are presented to show the accuracy of this hybrid method. Firstly, a straight beam with the following parameters is used in the calculation: beam length=
L
FEM
+
L
FSM
=
1
m , in which L
FEM and L
FSM represent the length of the FEM and FSM zones on the beam structure; Young’s modulus E
=
2.1
×
10
11
Pa; the second moment of inertia of the cross-section I
=
π
r
4
/
4 with r
=
0.01
m . The Fourier series are truncated to M
=
20 for the FSM, while the element number for the FEM is set as 50. In order to study the influence of the relative proportion of these two subdomains on the prediction results, the length of the FSM zone is varied from 0 to 1 with the increment of 0.2 while the overall length of the beam is kept unchanged as =
1
m . The results are given in Table 1 for a beam with clamped-clamped boundary condition, which can be readily achieved by setting all the end restraint stiffness into a very large number (the translational springs are k
w
0
=
5
×
1
0
7, k
w
1
=
5
×
1
0
7 while the rotation springs at the boundary are K
w
0
=
5
×
1
0
7, K
w
1
=
5
×
1
0
7 and the coupled springs are k
c
=
5
×
1
0
7, K
c
=
5
×
1
0
7) for FSM zone side and fixing the end node on FEM zone side. As can be observed, the relative proportion between these two computational zones almost has no effect on the accuracy of the predicted frequencies. Such good results indicate that the proposed hybrid method is effective.Table 1
Effect of the division ratio of the two computational domains on the modal frequencies of the beams structure with clamped-clamped boundary condition.
Mode order
μ
i
=
L
/
π
(
ω
i
ρ
A
/
EI
)
1
/
2, M = 20, E
n = 50
L
FSM (m)
Analytical results
1.0
0.8
0.6
0.4
0.2
0.0
L
FEM (m)
0.0
0.2
0.4
0.6
0.8
1.0
1
1.5056
1.5056
1.5056
1.5056
1.5056
1.5056
1.5055
2
2.4998
2.4998
2.4997
2.4998
2.4997
2.4998
2.4996
3
3.5000
3.5000
3.5000
3.5000
3.5000
3.5000
3.5000
4
4.5000
4.5000
4.5000
4.5000
4.5000
4.5000
4.4999
5
5.5000
5.5000
5.4999
5.4999
5.4999
5.4999
5.5000
6
6.5000
6.5000
6.4998
6.5000
6.5000
6.4998
6.5001
7
7.5000
7.5000
7.4999
7.4998
7.4999
7.4999
7.5002
8
8.5000
8.5000
8.5000
8.4996
8.4997
8.5001
8.5003
9
9.5000
9.5000
9.4999
9.5000
9.5000
9.5000
9.5004
10
10.5000
10.4999
10.4993
10.4992
10.4992
10.4995
10.5007Table2 compares the efficiency and accuracy of the FEM, the FSM, and the hybrid method. All these three methods are implemented on the same computational platform in MATLAB. A relative high frequency range up to the 20th mode is computed by these three approaches. Firstly, FEM and FSM are used to make prediction and it can be seen that FSM is much more efficient than FEM, since its dimension of final system matrix is much smaller than that of FEM. For the proposed hybrid method, it can be observed again that the prediction is accurate enough for capturing the modal characteristics of this free-free beam.Table 2
Comparison of computing time of the hybrid method with pure FEM and FEM with different truncated Fourier terms and/or FEM element numbers, in which a completely free-free beam structure is used.
Method
FSM truncatedterm, M
FEM elementnumber, E
n
The 20th modal frequency
Degrees of freedom
Computing time (s)
Analytical
—
—
19.50000000
—
—
FSM
80
—
19.50000044
85
0.016
FEM
—
600
19.50000073
1202
9.235
Hybrid
40
300
19.49999961
647
7.66
FSM
60
—
19.5000032
65
0.016
FEM
—
500
19.50000153
1002
5.133
Hybrid
30
250
19.50000005
537
3.947
FSM
50
—
19.5000114
55
0.015
FEM
—
250
19.50002436
502
0.64
Hybrid
25
125
19.50001615
282
0.499
FSM
30
—
19.50046727
35
<0.01
FEM
—
100
19.50094026
202
0.047
Hybrid
15
50
19.50069514
122
0.031In the first three lines of Table2, we see that in these three cases the accuracy class of the 20th modal frequency calculated is equal. The degrees of freedom of the matrix used for calculation are 85, 600, and 345 for FSM, FEM, and current hybrid method, respectively. Although the CPU times are not a good measure of an algorithm’s performance as they do not provide much insight about an algorithm’s behaviour and they are still presented here as references which are 0.016(s), 9.235(s), and 7.66(s), respectively, for the first three cases. Compared with FEM, since the FSM with rapid convergence and high efficiency is incorporated into this structure, the overall degrees of freedom will be reduced considerably, say, from 600 to 345 for the first three cases. The CPU time consumed is also reduced compared with pure FEM. The remaining cases in this table are all the same: the current method would reduce DOF significantly compared with pure FEM while the accuracy class of the frequency calculated is still kept in the same level.


## 3. Theoretical Formulation of 2-D Structures

### 3.1. Governing Equation of Rectangular Plates

Similar to the framework for vibration analysis of beam structures, vibrating plate structures with elastic edge supports can also be dealt with by the energy principle. For a rectangular plate with widtha and length b, the total potential energy V and total kinetic energy T are given by
(11)
V
=
D
2
∫
0
a
∫
0
b
(
(
∂
2
w
∂
x
2
)
2
+
(
∂
2
w
∂
y
2
)
2
)
d
x
d
y
+
D
2
∫
0
a
∫
0
b
(
2
ν
∂
2
w
∂
x
2
∂
2
w
/
∂
y
2
G
G
G
G
G
h
G
G
G
G
+
2
(
1
-
ν
)
(
∂
2
w
∂
x
∂
y
)
2
)
d
x
d
y
+
1
2
∫
0
b
(
k
x
0
w
2
+
K
x
0
(
∂
w
∂
x
)
2
)
x
=
0
d
y
+
1
2
∫
0
b
(
k
x
1
w
2
+
K
x
1
(
∂
w
∂
x
)
2
)
x
=
a
d
y
+
1
2
∫
0
a
(
k
y
0
w
2
+
K
y
0
(
∂
w
∂
y
)
2
)
y
=
0
d
x
+
1
2
∫
0
a
(
k
y
1
w
2
+
K
y
1
(
∂
w
∂
y
)
2
)
y
=
b
d
x
,
T
=
1
2
∫
0
a
∫
0
b
ρ
h
(
∂
w
∂
t
)
2
d
x
d
y
,
where w is the flexural displacement field, and D, ν, ρ, and h are, respectively, the flexural rigidity, the Poisson ratio, the mass density, and the plate thickness. k
x
0 and k
x
a (k
y
0 and k
y
b) are the constants of linearly distributed springs and K
x
0 and K
x
a (K
y
0 and K
y
b) the rotational spring constants at x
=
0 and x
=
a (y
=
0 and y
=
b), respectively.In the FEM description of a plate structure, each node has a lateral displacement and two rotations. The equation for vibration analysis of a plate structure can also be obtained in the form of (2).In this study, the displacement function of a rectangular plate for the FSM will be sought in the form of series expansions as follows [24]:
(12)
w
(
x
,
y
)
=
∑
m
=
0
∞
∑
n
=
0
∞
A
m
n
cos
⁡
λ
a
m
x
cos
⁡
λ
b
n
y
+
∑
l
=
1
4
(
ξ
l
b
(
y
)
∑
m
=
0
∞
c
m
l
cos
⁡
λ
a
m
x
h
h
h
h
h
h
h
h
+
ξ
l
a
(
x
)
∑
n
=
0
∞
d
n
l
cos
⁡
λ
b
n
y
)
,
where λ
a
m
=
m
π
/
a, λ
b
n
=
n
π
/
b, and ξ
j
a (or ξ
j
b) represent a set of closed-form supplementary functions, defined as
(13)
ξ
1
a
(
x
)
=
9
a
4
π
sin
(
π
x
2
a
)
-
a
12
π
sin
(
3
π
x
2
a
)
,
(14)
ξ
2
a
(
x
)
=
-
9
a
4
π
cos
⁡
(
π
x
2
a
)
-
a
12
π
cos
⁡
(
3
π
x
2
a
)
,
(15)
ξ
3
a
(
x
)
=
a
3
π
3
sin
(
π
x
2
a
)
-
a
3
3
π
3
sin
(
3
π
x
2
a
)
,
(16)
ξ
4
a
(
x
)
=
-
a
3
π
3
cos
⁡
(
π
x
2
a
)
-
a
3
3
π
3
cos
⁡
(
3
π
x
2
a
)
.
The form of ξ
j
b is much the same as ξ
j
a except that a should be replaced by b while x should be replaced by y. A procedure similar to that for the beam in Section 2 can be established to determine all the coefficients in (12), in which the series is truncated to m
=
M and n
=
N.


### 3.2. Combination of the FEM and the FSM

Compared with the hybrid beam model, the coefficients which couple these two methods are a little different. As illustrated in Figure2, a plate structure is divided into two zones: the left is handled by the FEM and the right by the FSM. Geometrically, in the coupled beam model, the two different beam zones are connected point to point, while for the coupled plate the two zones are connected line to line. As the rectangular plate is one type of structure which possesses symmetric characteristics, the edge on the left side of the FSM plate is considered here. For this paper, the model for FSM is based on classical Poisson-Kirchhoff theory, which requires C1 continuity in the framework of the finite element method [1]. Certain types of elements based on Poisson-Kirchhoff theory have compatibility problems as the rotational displacements would become discontinuous between adjacent elements on the edge. Hence, for the sake of generality, the transverse springs between these two zones are set to be linearly distributed while the rotation springs are concentrated to nodes. For the FSM zone, the displacement and rotation of the edge at x
=
0 can be expressed as
(17)
w
(
x
,
y
)
|
x
=
0
=
∑
m
=
0
∞
∑
n
=
0
∞
A
m
n
cos
⁡
λ
b
n
y
+
∑
l
=
1
4
(
ξ
l
b
(
y
)
∑
m
=
0
∞
c
m
l
j
j
j
j
j
j
j
j
j
j
j
j
+
ξ
l
a
(
0
)
∑
n
=
0
∞
d
n
l
cos
⁡
λ
b
n
y
)
,
∂
w
(
x
,
y
)
∂
x
|
x
=
0
=
∑
n
=
0
∞
d
n
1
cos
⁡
λ
b
n
y
.Figure 2
Schematic of elastically restrained plate structure divided into FEM zone and FSM zone, with two types of fictitious springs introduced at zone interface.For then
e element bordering the interface, the four nodes of the element are defined as i
n
e
,
j
n
e
,
l
n
e
,
and
m
n
e.And then the potential energy stored in linear spring can now be calculated as(18)
V
l
c
=
1
2
∫
y
n
e
j
y
n
e
m
k
l
c
(
w
n
e
|
x
=
0
-
w
f
s
|
x
=
0
)
2
d
y
,
where y
n
e
j and y
n
e
m are the coordinates of the nodes j
n
eand m
n
e, respectively, which are located on the interface line.The type of element used here is a 4-node incompatible plate element, in which the normal derivative at the edge is discontinuous. As a result, the rotation displacement at the edge cannot be obtained directly from derivatives of (12). For overcoming difficulty, one procedure is introduced here taking care of the uniformly distributed rotational spring along the interface line: they are replaced by a certain number of discrete point-to-point rotational springs between each pair of the nodes having the same spatial coordinates on either side of the interface line. Suppose that the number of the pairs of such nodes is Q
c
p. For these discrete rotation springs, the potential energy stored in them can be expressed as
(19)
V
r
c
=
1
2
∑
n
c
p
=
1
Q
c
p
K
r
c
(
θ
x
n
c
p
-
θ
f
s
|
x
=
0
,
y
=
y
n
c
p
)
2
,
where y
n
c
p represents the y-coordinate of the n
c
p nodes on the interface line and θ
x
n
c
p is the rotations about the x-axis at the n
c
p nodes.Substituting (17) into (18) and (19) and then setting the derivatives of the resultant equations (with respect to each unknown Fourier coefficient and unknown displacement w
n and θ
n) to zero, the coupled matrix K
couple can be obtained. In this way, the overall system equation can be written as
(20)
(
[
K
FEM
0
0
K
FSM
]
+
[
K
couple
]
JJJJ
-
ω
2
[
M
FEM
0
0
M
FSM
]
)
{
Λ
Ψ
}
=
{
0
}
,
in which
(21)
Λ
=
{
w
1
,
θ
x
1
,
θ
y
1
,
w
2
,
θ
x
2
,
θ
y
2
,
…
,
w
Q
,
θ
x
Q
,
θ
y
Q
}
T
,
Ψ
=
{
A
00
,
A
01
…
A
M
1
…
,
A
1
N
…
A
M
N
,
H
H
c
11
…
c
M
2
,
…
c
M
4
,
d
11
…
d
N
2
,
…
d
N
4
}
T
.
As the detailed coupled matrices can also be easily obtained according to the same procedures the beam coupled matrices are obtained and details of (20) are not given here.


### 3.3. Numerical Results and Discussion

In this section, the accuracy and effectiveness of the proposed hybrid method will be demonstrated for predicting the modal behaviour of a 2-D plate structure with various boundary conditions and/or different shapes. As the first example, a rectangular plate with simply supported boundary condition at all its edges is considered. For the FSM, a simply supported edge can be treated as a special case while the stiffness coefficients for the translational springs become infinitely large and the rotational spring becomes zero. The plate dimension is 5 m wide and 10 m long with thicknessh
=
0
.
008
m; the material parameters used for the plate are E
=
2.1
×
10
11
Pa and ρ
=
7800
Kg
/
m
3. As mentioned above, in numerical calculations the series expansion in (13) has to be truncated to M
=
N
=
12. The plate is divided by a straight line at the centre of the length into two square zones with all the width and length equal to 5 m. For the FEM model, 20
×
40 elements are chosen so that higher (such as the 8th, 9th, and 10th) frequencies given in Table 3 would then become accurate enough for comparison. Table 3 shows the first ten natural frequencies of a rectangular plate simply supported at all its edges (the translational springs are k
x
0
=
0
,
k
x
a
,
k
y
0
,
k
y
b
=
5
×
1
0
7 while the rotation springs at the boundary are K
x
0
=
0, K
x
a
,
K
y
0, K
y
b
=
5
×
1
0
7 and the coupled springs k
l
c
=
5
×
1
0
7
,
K
r
c
=
5
×
1
0
7). Several mode shapes are also given in Figure 3.Table 3
Comparison of the first ten natural frequencies of the rectangular plate.
Mode order
Analytical
FEAa
Hybridb
Difference
DOF = 2583
DOF = 1357
—
1
0.9866
0.98624
0.9868
−0.06%
2
1.5785
1.5774
1.5796
−0.14%
3
2.5651
2.5624
2.5645
−0.08%
4
3.3544
3.3522
3.354
−0.05%
5
3.9463
3.9411
3.9447
−0.09%
6
3.9463
3.9411
3.9458
−0.12%
7
4.9329
4.9223
4.9275
−0.11%
8
5.7221
5.7131
5.7177
−0.08%
9
6.3141
6.2958
6.3057
−0.16%
10
7.3007
7.2930
7.2989
−0.08%
Note:aFEA results are calculated by using ANSYS with 20 × 40 elements.
bResults with M = N = 12 at the FSM zone and 20 × 20 elements at the FEM zone.Plots of the mode shape (meshed line denotes the FEM zone, while the continuous part denotes the FSM zone) for the rectangular. The (a) first, (b) second, (c) third, (d) fourth, (e) sixth, and (f) tenth mode shapes.
(a)
(b)
(c)
(d)
(e)
(f)It is possible that the matrixK becomes ill-conditioned if the coupled springs, k
l
c and K
r
c, are set to very large numbers. In Figure 4, the first 6 mode frequencies of the rectangular plate are given with the coupled springs, k
l
c and K
r
c, varied from 0 to 1010. It can be clearly seen from this figure that the influence of the coupled springs on the frequencies calculated is minimal when the springs changed within a wide range from 106 to 1010.Figure 4
Influence of coupled springs on the former 6 mode frequencies of the rectangular plate.In order to verify the accuracy and stability of the proposed hybrid method thoroughly, frequencies of higher modes are provided in Table4. It can be seen clearly from this table that the hybrid method is slightly more accurate than the FEM compared with analytical results. However, because the FEM is known to be incapable of producing good results at very high frequencies and the hybrid method uses the FEM, it is expected that the hybrid method will produce more accurate results of high frequencies than the FEM but still cannot produce very accurate results of very high frequencies.Table 4
Comparison of higher natural frequencies of the rectangular plate.
Mode order
Analytical
FEAa
ERROR
Hybridb
ERROR
10
7.3007
7.2930
0.10%
7.2989
0.02%
20
12.8255
12.8055
0.16%
12.8225
0.02%
30
19.1396
19.0353
0.54%
19.0681
0.37%
40
23.0859
22.8829
0.88%
22.9219
0.71%
50
29.4000
29.0709
1.12%
29.1143
0.97%
60
35.5168
35.1455
1.05%
35.1161
1.13%
70
40.4496
40.0103
1.09%
40.1129
0.83%
80
45.7772
45.2644
1.12%
45.3688
0.89%
90
51.3020
51.0894
0.41%
51.2793
0.04%
100
56.8268
55.60429
2.15%
55.7703
1.86%
200
109.9047
106.3049
3.28%
106.4953
3.10%
Note:aFEA results are calculated by using ANSYS with 20 × 40 elements.
bResults with M = N = 12 at the FSM zone and 20 × 20 elements at the FEM zone.Another example of a rectangular plate (identical to the previous FSM rectangular plate with its width and length equal to 5 m) plus a matching semicircular plate (with its diameter equal to the width of the rectangular plate) is also solved and the results are presented in Table5 with their corresponding mode shapes plotted in Figure 5. For the semicircular plate, the finite element matrices for bending are directly extracted from the ANSYS software based on the shell 63 element. As shown in Figure 5, elements used in this example are irregular- or even triangular-shaped. And the results show the generality of the present method when using finite elements with irregular shapes.Table 5
Comparison of the first ten natural frequencies of the rectangular-half circular plate.
Mode order
FEAa
Hybridb
Difference
DOF = 2472
DOF = 781
—
1
0.8266
0.8265
0.01%
2
0.8309
0.8311
−0.02%
3
1.876
1.876
0.00%
4
1.9316
1.9308
0.04%
5
2.1311
2.1299
0.06%
6
2.8449
2.8444
0.02%
7
3.5087
3.5053
0.10%
8
3.9185
3.9167
0.05%
9
4.7067
4.7043
0.05%
10
4.9564
4.9553
0.02%
Note:aFEA results are calculated by using ANSYS with 773 elements.
bResults with M = N = 12 at the FSM zone and 22 elements at the edge of the coupled line and the element matrix of the semicircular plate were derived from the ANSYS software.Plots of the mode shape (meshed line denotes the FEM zone, while the continuous part denotes the FSM zone) for the rectangular-circular plate. The (a) first, (b) second, (c) third, (d) fourth, (e) sixth, and (f) twelfth mode shapes.
(a)
(b)
(c)
(d)
(e)
(f)It can also be seen from this example that the present method can be applied to structures with geometrically irregular boundary conditions. In these two 2-D examples, the validity and accuracy of the presented hybrid method are also confirmed by comparison with FEM results.


## 3.1. Governing Equation of Rectangular Plates

Similar to the framework for vibration analysis of beam structures, vibrating plate structures with elastic edge supports can also be dealt with by the energy principle. For a rectangular plate with widtha and length b, the total potential energy V and total kinetic energy T are given by
(11)
V
=
D
2
∫
0
a
∫
0
b
(
(
∂
2
w
∂
x
2
)
2
+
(
∂
2
w
∂
y
2
)
2
)
d
x
d
y
+
D
2
∫
0
a
∫
0
b
(
2
ν
∂
2
w
∂
x
2
∂
2
w
/
∂
y
2
G
G
G
G
G
h
G
G
G
G
+
2
(
1
-
ν
)
(
∂
2
w
∂
x
∂
y
)
2
)
d
x
d
y
+
1
2
∫
0
b
(
k
x
0
w
2
+
K
x
0
(
∂
w
∂
x
)
2
)
x
=
0
d
y
+
1
2
∫
0
b
(
k
x
1
w
2
+
K
x
1
(
∂
w
∂
x
)
2
)
x
=
a
d
y
+
1
2
∫
0
a
(
k
y
0
w
2
+
K
y
0
(
∂
w
∂
y
)
2
)
y
=
0
d
x
+
1
2
∫
0
a
(
k
y
1
w
2
+
K
y
1
(
∂
w
∂
y
)
2
)
y
=
b
d
x
,
T
=
1
2
∫
0
a
∫
0
b
ρ
h
(
∂
w
∂
t
)
2
d
x
d
y
,
where w is the flexural displacement field, and D, ν, ρ, and h are, respectively, the flexural rigidity, the Poisson ratio, the mass density, and the plate thickness. k
x
0 and k
x
a (k
y
0 and k
y
b) are the constants of linearly distributed springs and K
x
0 and K
x
a (K
y
0 and K
y
b) the rotational spring constants at x
=
0 and x
=
a (y
=
0 and y
=
b), respectively.In the FEM description of a plate structure, each node has a lateral displacement and two rotations. The equation for vibration analysis of a plate structure can also be obtained in the form of (2).In this study, the displacement function of a rectangular plate for the FSM will be sought in the form of series expansions as follows [24]:
(12)
w
(
x
,
y
)
=
∑
m
=
0
∞
∑
n
=
0
∞
A
m
n
cos
⁡
λ
a
m
x
cos
⁡
λ
b
n
y
+
∑
l
=
1
4
(
ξ
l
b
(
y
)
∑
m
=
0
∞
c
m
l
cos
⁡
λ
a
m
x
h
h
h
h
h
h
h
h
+
ξ
l
a
(
x
)
∑
n
=
0
∞
d
n
l
cos
⁡
λ
b
n
y
)
,
where λ
a
m
=
m
π
/
a, λ
b
n
=
n
π
/
b, and ξ
j
a (or ξ
j
b) represent a set of closed-form supplementary functions, defined as
(13)
ξ
1
a
(
x
)
=
9
a
4
π
sin
(
π
x
2
a
)
-
a
12
π
sin
(
3
π
x
2
a
)
,
(14)
ξ
2
a
(
x
)
=
-
9
a
4
π
cos
⁡
(
π
x
2
a
)
-
a
12
π
cos
⁡
(
3
π
x
2
a
)
,
(15)
ξ
3
a
(
x
)
=
a
3
π
3
sin
(
π
x
2
a
)
-
a
3
3
π
3
sin
(
3
π
x
2
a
)
,
(16)
ξ
4
a
(
x
)
=
-
a
3
π
3
cos
⁡
(
π
x
2
a
)
-
a
3
3
π
3
cos
⁡
(
3
π
x
2
a
)
.
The form of ξ
j
b is much the same as ξ
j
a except that a should be replaced by b while x should be replaced by y. A procedure similar to that for the beam in Section 2 can be established to determine all the coefficients in (12), in which the series is truncated to m
=
M and n
=
N.


## 3.2. Combination of the FEM and the FSM

Compared with the hybrid beam model, the coefficients which couple these two methods are a little different. As illustrated in Figure2, a plate structure is divided into two zones: the left is handled by the FEM and the right by the FSM. Geometrically, in the coupled beam model, the two different beam zones are connected point to point, while for the coupled plate the two zones are connected line to line. As the rectangular plate is one type of structure which possesses symmetric characteristics, the edge on the left side of the FSM plate is considered here. For this paper, the model for FSM is based on classical Poisson-Kirchhoff theory, which requires C1 continuity in the framework of the finite element method [1]. Certain types of elements based on Poisson-Kirchhoff theory have compatibility problems as the rotational displacements would become discontinuous between adjacent elements on the edge. Hence, for the sake of generality, the transverse springs between these two zones are set to be linearly distributed while the rotation springs are concentrated to nodes. For the FSM zone, the displacement and rotation of the edge at x
=
0 can be expressed as
(17)
w
(
x
,
y
)
|
x
=
0
=
∑
m
=
0
∞
∑
n
=
0
∞
A
m
n
cos
⁡
λ
b
n
y
+
∑
l
=
1
4
(
ξ
l
b
(
y
)
∑
m
=
0
∞
c
m
l
j
j
j
j
j
j
j
j
j
j
j
j
+
ξ
l
a
(
0
)
∑
n
=
0
∞
d
n
l
cos
⁡
λ
b
n
y
)
,
∂
w
(
x
,
y
)
∂
x
|
x
=
0
=
∑
n
=
0
∞
d
n
1
cos
⁡
λ
b
n
y
.Figure 2
Schematic of elastically restrained plate structure divided into FEM zone and FSM zone, with two types of fictitious springs introduced at zone interface.For then
e element bordering the interface, the four nodes of the element are defined as i
n
e
,
j
n
e
,
l
n
e
,
and
m
n
e.And then the potential energy stored in linear spring can now be calculated as(18)
V
l
c
=
1
2
∫
y
n
e
j
y
n
e
m
k
l
c
(
w
n
e
|
x
=
0
-
w
f
s
|
x
=
0
)
2
d
y
,
where y
n
e
j and y
n
e
m are the coordinates of the nodes j
n
eand m
n
e, respectively, which are located on the interface line.The type of element used here is a 4-node incompatible plate element, in which the normal derivative at the edge is discontinuous. As a result, the rotation displacement at the edge cannot be obtained directly from derivatives of (12). For overcoming difficulty, one procedure is introduced here taking care of the uniformly distributed rotational spring along the interface line: they are replaced by a certain number of discrete point-to-point rotational springs between each pair of the nodes having the same spatial coordinates on either side of the interface line. Suppose that the number of the pairs of such nodes is Q
c
p. For these discrete rotation springs, the potential energy stored in them can be expressed as
(19)
V
r
c
=
1
2
∑
n
c
p
=
1
Q
c
p
K
r
c
(
θ
x
n
c
p
-
θ
f
s
|
x
=
0
,
y
=
y
n
c
p
)
2
,
where y
n
c
p represents the y-coordinate of the n
c
p nodes on the interface line and θ
x
n
c
p is the rotations about the x-axis at the n
c
p nodes.Substituting (17) into (18) and (19) and then setting the derivatives of the resultant equations (with respect to each unknown Fourier coefficient and unknown displacement w
n and θ
n) to zero, the coupled matrix K
couple can be obtained. In this way, the overall system equation can be written as
(20)
(
[
K
FEM
0
0
K
FSM
]
+
[
K
couple
]
JJJJ
-
ω
2
[
M
FEM
0
0
M
FSM
]
)
{
Λ
Ψ
}
=
{
0
}
,
in which
(21)
Λ
=
{
w
1
,
θ
x
1
,
θ
y
1
,
w
2
,
θ
x
2
,
θ
y
2
,
…
,
w
Q
,
θ
x
Q
,
θ
y
Q
}
T
,
Ψ
=
{
A
00
,
A
01
…
A
M
1
…
,
A
1
N
…
A
M
N
,
H
H
c
11
…
c
M
2
,
…
c
M
4
,
d
11
…
d
N
2
,
…
d
N
4
}
T
.
As the detailed coupled matrices can also be easily obtained according to the same procedures the beam coupled matrices are obtained and details of (20) are not given here.


## 3.3. Numerical Results and Discussion

In this section, the accuracy and effectiveness of the proposed hybrid method will be demonstrated for predicting the modal behaviour of a 2-D plate structure with various boundary conditions and/or different shapes. As the first example, a rectangular plate with simply supported boundary condition at all its edges is considered. For the FSM, a simply supported edge can be treated as a special case while the stiffness coefficients for the translational springs become infinitely large and the rotational spring becomes zero. The plate dimension is 5 m wide and 10 m long with thicknessh
=
0
.
008
m; the material parameters used for the plate are E
=
2.1
×
10
11
Pa and ρ
=
7800
Kg
/
m
3. As mentioned above, in numerical calculations the series expansion in (13) has to be truncated to M
=
N
=
12. The plate is divided by a straight line at the centre of the length into two square zones with all the width and length equal to 5 m. For the FEM model, 20
×
40 elements are chosen so that higher (such as the 8th, 9th, and 10th) frequencies given in Table 3 would then become accurate enough for comparison. Table 3 shows the first ten natural frequencies of a rectangular plate simply supported at all its edges (the translational springs are k
x
0
=
0
,
k
x
a
,
k
y
0
,
k
y
b
=
5
×
1
0
7 while the rotation springs at the boundary are K
x
0
=
0, K
x
a
,
K
y
0, K
y
b
=
5
×
1
0
7 and the coupled springs k
l
c
=
5
×
1
0
7
,
K
r
c
=
5
×
1
0
7). Several mode shapes are also given in Figure 3.Table 3
Comparison of the first ten natural frequencies of the rectangular plate.
Mode order
Analytical
FEAa
Hybridb
Difference
DOF = 2583
DOF = 1357
—
1
0.9866
0.98624
0.9868
−0.06%
2
1.5785
1.5774
1.5796
−0.14%
3
2.5651
2.5624
2.5645
−0.08%
4
3.3544
3.3522
3.354
−0.05%
5
3.9463
3.9411
3.9447
−0.09%
6
3.9463
3.9411
3.9458
−0.12%
7
4.9329
4.9223
4.9275
−0.11%
8
5.7221
5.7131
5.7177
−0.08%
9
6.3141
6.2958
6.3057
−0.16%
10
7.3007
7.2930
7.2989
−0.08%
Note:aFEA results are calculated by using ANSYS with 20 × 40 elements.
bResults with M = N = 12 at the FSM zone and 20 × 20 elements at the FEM zone.Plots of the mode shape (meshed line denotes the FEM zone, while the continuous part denotes the FSM zone) for the rectangular. The (a) first, (b) second, (c) third, (d) fourth, (e) sixth, and (f) tenth mode shapes.
(a)
(b)
(c)
(d)
(e)
(f)It is possible that the matrixK becomes ill-conditioned if the coupled springs, k
l
c and K
r
c, are set to very large numbers. In Figure 4, the first 6 mode frequencies of the rectangular plate are given with the coupled springs, k
l
c and K
r
c, varied from 0 to 1010. It can be clearly seen from this figure that the influence of the coupled springs on the frequencies calculated is minimal when the springs changed within a wide range from 106 to 1010.Figure 4
Influence of coupled springs on the former 6 mode frequencies of the rectangular plate.In order to verify the accuracy and stability of the proposed hybrid method thoroughly, frequencies of higher modes are provided in Table4. It can be seen clearly from this table that the hybrid method is slightly more accurate than the FEM compared with analytical results. However, because the FEM is known to be incapable of producing good results at very high frequencies and the hybrid method uses the FEM, it is expected that the hybrid method will produce more accurate results of high frequencies than the FEM but still cannot produce very accurate results of very high frequencies.Table 4
Comparison of higher natural frequencies of the rectangular plate.
Mode order
Analytical
FEAa
ERROR
Hybridb
ERROR
10
7.3007
7.2930
0.10%
7.2989
0.02%
20
12.8255
12.8055
0.16%
12.8225
0.02%
30
19.1396
19.0353
0.54%
19.0681
0.37%
40
23.0859
22.8829
0.88%
22.9219
0.71%
50
29.4000
29.0709
1.12%
29.1143
0.97%
60
35.5168
35.1455
1.05%
35.1161
1.13%
70
40.4496
40.0103
1.09%
40.1129
0.83%
80
45.7772
45.2644
1.12%
45.3688
0.89%
90
51.3020
51.0894
0.41%
51.2793
0.04%
100
56.8268
55.60429
2.15%
55.7703
1.86%
200
109.9047
106.3049
3.28%
106.4953
3.10%
Note:aFEA results are calculated by using ANSYS with 20 × 40 elements.
bResults with M = N = 12 at the FSM zone and 20 × 20 elements at the FEM zone.Another example of a rectangular plate (identical to the previous FSM rectangular plate with its width and length equal to 5 m) plus a matching semicircular plate (with its diameter equal to the width of the rectangular plate) is also solved and the results are presented in Table5 with their corresponding mode shapes plotted in Figure 5. For the semicircular plate, the finite element matrices for bending are directly extracted from the ANSYS software based on the shell 63 element. As shown in Figure 5, elements used in this example are irregular- or even triangular-shaped. And the results show the generality of the present method when using finite elements with irregular shapes.Table 5
Comparison of the first ten natural frequencies of the rectangular-half circular plate.
Mode order
FEAa
Hybridb
Difference
DOF = 2472
DOF = 781
—
1
0.8266
0.8265
0.01%
2
0.8309
0.8311
−0.02%
3
1.876
1.876
0.00%
4
1.9316
1.9308
0.04%
5
2.1311
2.1299
0.06%
6
2.8449
2.8444
0.02%
7
3.5087
3.5053
0.10%
8
3.9185
3.9167
0.05%
9
4.7067
4.7043
0.05%
10
4.9564
4.9553
0.02%
Note:aFEA results are calculated by using ANSYS with 773 elements.
bResults with M = N = 12 at the FSM zone and 22 elements at the edge of the coupled line and the element matrix of the semicircular plate were derived from the ANSYS software.Plots of the mode shape (meshed line denotes the FEM zone, while the continuous part denotes the FSM zone) for the rectangular-circular plate. The (a) first, (b) second, (c) third, (d) fourth, (e) sixth, and (f) twelfth mode shapes.
(a)
(b)
(c)
(d)
(e)
(f)It can also be seen from this example that the present method can be applied to structures with geometrically irregular boundary conditions. In these two 2-D examples, the validity and accuracy of the presented hybrid method are also confirmed by comparison with FEM results.


## 4. Conclusions

In this paper, a hybrid FEM-FSM modelling method for the vibration analysis of structures with elastic boundary conditions is proposed. The key point of this method is to divide the whole structure into FEM and FSM zones, and then fictitious springs are introduced to couple the two zones at their interface within the framework of energy principle. Comparative studies of the hybrid method with the FSM and FEM are carried out for modal characteristics, which confirm both the validity of the approach and its efficiency. The number of the degrees of freedom can also be significantly reduced compared with the conventional FEM. Although the concept of such a hybrid method is demonstrated only for a 1-D structure (Euler-Bernoulli beam) and two 2-D structures (plate) in this paper, the proposed approach can easily be extended to more complicated 3-D structures such as box-type structures, as well as structures with holes or notches.



---

*Source: 834834-2014-09-25.xml*