Calculation of wave propagation characteristics in pre-deformed periodic
lattice frame structures via Spectral Element Method
Marius MELLMANN1∗ ; Chuanzeng ZHANG1
1 Chair of Structural Mechanics, University of Siegen, Germany

Abstract
In recent years, the design, optimization and application of periodic materials and structures, such as phononic
and photonic crystals, are becoming popular in many engineering ﬁelds. For instance, the phononic crystals canmanipulate the elastic and acoustic wave propagation characteristics. Consequently, it is feasible to create theso-called band-gaps, which are certain frequency ranges, in which the propagation of elastic or acoustic waves
is prohibited (1, 2). Accordingly, such structures have many promising applications in sound and vibration
insulation as well as wave ﬁlters. In this paper, the dispersion curves of pre-deformed periodic lattice frame
structures are calculated using the Spectral Element Method (SEM), which is based on the exact solution of theelastic wave equations. Numerical examples will be presented and compared with FEM calculations.
Keywords: Periodic lattice structures, Phononic crystals, Spectral Element Method (SEM)

1

INTRODUCTION

For periodic lattice frame structures, it is possible to manipulate the vibration and elastic wave propagation
properties by altering the stiffness and density of the structural members or by adding local resonators. In
addition to the change of the stiffness and mass of the system,
it is also possible to control the vibration and wave propagation
properties by means of a speciﬁc pre-deformation of the strucPre-deformationtural elements, such as shown in Figure 1 (3, 4). The vibration
and wave propagation characteristics of such structures can beobtained by considering a representative unit-cell with imposed
Original latticeÄZig-zag³lattice
Bloch-Floquet boundary conditions and calculating the corresponding dispersion curves. This calculation can be carried out
Figure 1. Pre-deformation of a periodic latby different numerical methods such as the FEM, FDM, Plane
tice structure.
Wave Expansion (PWE) Method and Layered Multiple Scattering (LMS) Method (5). In this paper the dispersion curves of pre-deformed periodic lattice frame structures
are calculated using the Spectral Element Method (SEM), which is based on the exact solution of the elastic
wave equations. This results in the advantage, that it is sufﬁcient to represent each beam with just one spectral
element and the eigenvalue calculation still yields the exact eigenfrequencies.

2

SPECTRAL ELEMENT METHODThe Spectral Element Method used in this paper is based on the approach by Doyle (6) which uses theexact solution of the elastic wave equations to derive spectral element matrices. This method is not to be
mistaken for the Spectral Element Method proposed by Patera (7) which is based on the FEM with high ordershape functions. The approach by Doyle begins by transforming the governing partial differential equations of
motion from the time domain into the frequency domain via Fourier-Transformation. Consequently, the time
variable is replaced by the frequency variable and the time-domain partial differential equations are transformed∗ Corresponding author: mellmann@bau.uni-siegen.de3205

into frequency-domain ordinary differential equations. Subsequently, these equations are solved exactly and byusing the exact wave solutions frequency-dependent dynamic shape functions are derived (8). By means of
the dynamic shape functions the exact dynamic stiffness matrix, called the spectral element matrix, is ﬁnallyformulated similarly to the procedure commonly used in FEM (8).
In the following sections the derivation of the spectral element matrices for the rod and the Euler-Bernoulli
beam are presented. The formulation is based on the force-displacement relation method. For other methods ofspectral element formulation, like the variational method or the state-vector equation method, we refer to (8).
2.1 Derivation of the spectral element matrix for a rod
The partial differential equation of motion for free longitudinal vibration of a rod is given by
EAu − ρAü = 0,(1)

where u(x,t) is the longitudinal displacement, E is the Young’s modulus, A is the cross-sectional area and ρ isthe mass density. The solution of Eq. 1 is assumed asu(x,t) =1 N−1∑ Un (x, ωn )eiωnt .
N n=0

(2)

Substituting the approach of Eq. 2 into Eq. 1 yields the eigenvalue problem
EAU  + ω 2 ρAU = 0(3)

for each discrete frequency ωn . The general solution of Eq. 3 is assumed to beU(x) = ae−ik(ω)x .

(4)

Substituting Eq. 4 in Eq. 3 yields the dispersion relation
k2 − kL2 = 0(5)

k1 = −k2 = kL ,(6)

with the two real roots
where kL is the wave number of the longitudinal wave mode and is given byρA
kL = ω
.
(7)
EA
Finally, by inserting the roots in Eq. 4 the general solution for a ﬁnite rod element of length L is given by
U(x) = a1 e−ikL x + a2 eikL x = e(x, ω)a
e(x, ω) = e−ikL x eikL x ,

with

(8)

a = [a1 a2 ]T .

(9)

The values of a depend on the spectral nodal displacements of the ﬁnite rod element, pictured in Figure 2,
which can be related to the displacement ﬁeld via
d = [U1 U2 ]T = [U(0) U(L)]T .
U1
N1

U2

N ( x)L

N ( x) = EAu 'Figure 2. Sign convention for the rod element.

3206

(10)

N2

Evaluating Eq. 8 at these nodes yields
d = [e(0, ω) e(L, ω)]T a = HR (ω)a


where
HR (ω) =1

e+ikL L1

e−ikL L

(11)


.

(12)

The displacement ﬁeld in the ﬁnite rod element can now be represented in terms of the nodal DOFs by eliminating the constant vector a from Eq. 8 by using Eq. 11
U(x) = NR (x, ω)d(13)

NR (x, ω) = e(x, ω)H−1
R (ω) = [NR1 NR2 ],
NR1 (x, ω) = csc(kL L) sin(kL (L − x)),

(14)

where

NR2 (x, ω) = csc(kL L) sin(kL x).
The spectral components of the axial force are related to U(x) by
N(x) = EAU  (x).

(15)

The equilibrium condition between the spectral nodal axial forces and the inner forces deﬁned by the strengthof the materials yields(16)
f c (ω) = [N1 N2 ]T = [−N(0) + N(L)]T .
Substituting Eq. 13 and Eq. 15 into the right-hand side of Eq. 16 gives
SR (ω)d = f c (ω),
in which SR (ω) is the spectral element matrix for the ﬁnite rod element and is given by
EA SR11 SR12
SR (ω) == STR (ω),
L SR21 SR22
where

SR11 = SR22 = (kL L) cot(kL L),
SR12 = SR21 = −(kL L) csc(kL L).

(17)

(18)

(19)

2.2 Derivation of the spectral element matrix for a Euler-Bernoulli beam
The partial differential equation of motion for free bending vibration of a Euler-Bernoulli beam is given by
EIw + ρAẅ = 0,(20)

where w(x,t) is the transverse displacement, E is the Young’s modulus, A is the cross-sectional area, I is thearea moment of inertia and ρ is the mass density. The solution of Eq. 20 is assumed asw(x,t) =1 N−1∑ Wn (x, ωn )eiωnt .
N n=0

(21)

Substituting the approach of Eq. 21 into equation Eq. 20 yields the eigenvalue problem
EIW  − ω 2 ρAW = 0

(22)

3207

for each discrete frequency ωn . The general solution of Eq. 22 is assumed to be
W (x) = ae−ik(ω)x .

(23)

Substituting Eq. 23 in Eq. 22 yields the dispersion relation
k4 − kF4 = 0(24)

with two real roots and two imaginary roots
k1 = −k2 = kF ,k3 = −k4 = ikF ,(25)

where kF is the wave number of the bending wave mode and is given by
kF =

√


ω

ρA
EI

1/4
.

(26)

Finally, by inserting the roots in Eq. 23 the general solution for a ﬁnite Euler-Bernoulli beam element of length
L is given by(27)
W (x) = a1 e−ikF x + a2 e−kF x + a3 eikF x + a4 ekF x = e(x, ω)a
with



e(x, ω) = e−ikF x e−kF x e+ikF x e+kF x ,a = [a1 a2 a3 a4 ]T .

(28)

The values of a depend on the spectral nodal displacements of the ﬁnite Euler-Bernoulli beam element, shownin Figure 3, which can be related to the displacement ﬁeld via
d = [W1 θ1 W2 θ2 ]T = [W (0) θ (0) W (L) θ (L)]T .
W2

W1

θ2

θ1
M1

(29)

Q1

Q( x)

Q2

Q( x)

M ( x)M2

L
M ( x) = EIw ''Q ( x) = − EIw '''Figure 3. Sign convention for the rod element.
Evaluating Eq. 27 at these nodes yields
d = e(0, ω) e (0, ω) e(L, ω) e (L, ω)
where

⎡

1
⎢ −ikF
HB (ω) = ⎢⎣ e−ikF L−ikF e−ikF L

1−kF
e−kF L−kF e−kF L

Ta = HB (ω)a,1
ikF

e+ikF LikF e+ikF L

1
kF

e+kF L

(30)
⎤
⎥
⎥.
⎦

(31)

kF e+kF L

By eliminating the constant vector a from Eq. 27 by using Eq. 30, we can represent the displacement ﬁeld in
the ﬁnite Euler-Bernoulli beam element in terms of the nodal DOFs with
W (x) = NB (x, ω)d

3208

(32)

whereNB (x, ω) = e(x, ω)H−1
B (ω) = [NB1 NB2 NB3 NB4 ],
NB1 (x, ω) = η −1 kF cos x − cos(L − x) cosh L − cos L cosh(L − x)
+ cosh x + sin(L − x) sinh L − sin L sinh(L − x) ,NB2 (x, ω) = η−1NB3 (x, ω) = η−1− cosh(L − x) sin L + cosh L sin(L − x) + sin x− cos(L − x) sinh L + cos L sinh(L − x) + sinh x ,kF cos(L − x) − cos x cosh L − cos L cosh x(33)+ cosh(L − x) + sin x sinh L − sin L sinh x ,NB4 (x, ω) = −η−1− cosh x sin L + cosh L sin x + sin(L − x)

+ cos x sinh L + cos L sinh x + sinh(L − x) ,
η = 2kF (1 − cos L cosh L),
x = kF x, L = kF L.
The spectral components of the bending moment and transverse shear force are related to W (x) by
Q(x) = −EIW  (x),M(x) = EIW  (x).

(34)

The equilibrium condition of the spectral nodal transverse shear forces and bending moments with the corresponding forces and moments deﬁned by the strength of the materials yields
f c (ω) = [Q1 M1 Q2 M2 ]T = [−Q(0) − M(0) + Q(L) + M(L)]T .

(35)

Substituting Eq. 32 and Eq. 34 into the right-hand side of Eq. 35 gives
SB (ω)d = f c (ω),(36)

in which SB (ω) is the spectral element matrix for the ﬁnite Euler-Bernoulli beam element and is given by⎡⎤
SB11 SB12 SB13 SB14
EI ⎢SB12 SB22 SB23 SB24 ⎥⎥ = STB (ω),
SB (ω) = 3 ⎢(37)
L ⎣SB13 SB23 SB33 SB34 ⎦SB14 SB24 SB34 SB44
where

3

SB11 = SB33 = ΔB L (cos L sinh L + sin L cosh L),
SB22 = SB44 = ΔB L kF−2 (− cos L sinh L + sin L cosh L),
3

SB12 = −SB34 = ΔB L kF−1 sin L sinh L,
3

3

SB13 = −ΔB L (sin L + sinh L),
SB14 = −SB23 = ΔB L kF−1 (− cos L + cosh L),
3(38)

SB24 = ΔB L kF−2 (− sin L + sinh L),
1,
ΔB =1 − cos L cosh L
L = kF L.
3

The extended spectral element matrix for the Euler-Bernoulli beam element including longitudinal vibrations canbe obtained by combining the spectral element matrices for the rod and the Euler-Bernoulli beam, yielding a
6 × 6 element matrix, which is used for the following calculations.

3209

3

PROBLEM FORMULATION FOR A UNIT-CELL

For an efﬁcient computation using the Spectral Element Method (SEM), an inﬁnitely large structure can bemapped by a representative unit-cell illustrated in Figure 4. The unit-cell of the undeformed lattice is depicted
in Figure 4a and the unit-cell of the pre-deformed "zig-zag" lattice is shown in Figure 4b. The periodicity of thestructure is accounted for by using the Bloch-Floquet theory. According to the Bloch-Floquet theory, so-called
H
H
c.)
a.)
b.)

H

Irreducible Brillouin zone
Μ

:

H

H

l1

l2

Γ

Χ

H

l1

Figure 4. Representative unit-cells of the a.) undeformed and b.) the pre-deformed, "zig-zag" shaped, periodic
lattice frame structure. c.) Corresponding irreducible Brillouin zone.
input and output nodes are deﬁned for the unit-cell. A simple 1-D example consisting of 2 spectral elements
is illustrated in Figure 5. The ﬁrst input node is deﬁned on the left side of the unit-cell with the DOFs
D1 , D2 , D3 and the corresponding output node is deﬁned on the right side of the unit-cell with the DOFs
D7 , D8 , D9 , respectively. Imposing Bloch-Floquet boundary conditions yields
Di,out = [D7 D8 D9 ]T = Di,in eik = [D1 D2 D3 ]T eik(39)

and accounts for the periodicity on the right side of the unit-cell. Respectively, the periodicity on the left sideof the unit-cell is considered by
Di,out = [D10 D11 D12 ]T = Di,in e−ik = [D4 D5 D6 ]T e−ik ,(40)

in which k = [kx e1 ky e2 ]T is the wave vector describing the phase shift of a wave between a unit-cell and it’sneighboring unit-cells. The wave vector, in turn, consists of the wave numbers in the x− and y−directions, each
multiplied by the associated lattice-tiling vectors e1 and e2 . In the 1-D case considered above the lattice-tiling
vector in y−direction vanishes and the wave vector is reduced to a scalar.
D3

D12
D11

D10

Master node
Slave node

D2

D9

D6
D1

D5

A1, I1, E1, ȡ1, l1

D4

D8

D7

A2, I2, E2, ȡ2, l2

Unit-cell

Figure 5. Application of Bloch-Floquet boundary conditions to 1-D unit-cell.
Imposing the boundary conditions reduces the size of the systems spectral matrix from 12 × 12 to 6 × 6. Finally,
the dispersion curves of the unit-cell, which describe the vibration and wave propagation characteristics of theoverall structure, can be determined. For this purpose the eigenvalue problem can be written as
Sges (ω, k)D = 0,

3210

(41)

which leads to a transcendental equation. This equation has to be solved for for a certain number of wave
vectors. The relevant wave vectors are deﬁned by the irreducible Brillouin zone in the reciprocal lattice, as
illustrated in Figure 4c. For a more detailed explanation of Bloch-Floquet theory, the Brillouin zone and thedetermination of dispersion curves the authors refer to (2, 9, 10).

4

EXAMPLE AND DISCUSSION

For the following example, a unit-cell with the sizes l1 = l2 = 100 mm and the cross-sectional values
A = 15 mm2 and Iy = 10 mm4 is analyzed. Furthermore, the material parameters are set to E = 30000 MPa
and ρ = 3000 kg/m3 . The ﬁrst 8 eigenfrequencies are calculated as functions of the wave vector yielding thecorresponding dispersion curves. Without pre-deformation, the calculation results in dispersion curves with noband-gaps in the considered frequency range, as shown in Figure 6a. Accordingly, this means that the vibrations
a.)

b.)

c.)

Band-gap: 4321-5879 Hz

Band-gap: 4338-5898 Hz

Band-gap: 2087-2939 Hz

Band-gap: 2087-2943 Hz

Figure 6. Dispersion curves of the a.) undeformed lattice, b.) the "zig-zag" lattice calculated by SEM (W =8 mm) and c.) the "zig-zag" lattice calculated by FEM with 10 ﬁnite elements per beam (W = 8 mm).
and elastic waves in the considered frequency range can propagate through the structure. In contrast, in caseof a "zig-zag" shaped pre-deformation of the structure with an amplitude W = 8 mm, the dispersion curves asshown in Figure 6b exhibit a gap between the 4th and 5th band as well as a gap between the 6th and 7th
band. Consequently, the vibrations and elastic waves with a frequency between 2087 and 2939 Hz, respectively
between 4321 and 5879 Hz, can notpropagate in the structure. Depending on the amplitude and the shapeof the pre-deformation, the dispersion curves as presented in Figure
6b can be further optimized, e.g.
wider or further band-gaps can beachieved, as shown in (11). The stationary solution for a lattice consisting of 20 × 20 unit-cells, excited by
the frequency f = 4600 Hz, is shownin Figure 7. It is obvious, that the
propagation of the wave is strongly
attenuated.
The dispersion curves for the sameunit-cell but calculated with FEM, in Figure 7. Stationary solution (by COMSOL) of the displacement ﬁeld
which each beam is discretized into for a pre-deformed lattice, with an amplitude of W = 8 mm, consisting
10 ﬁnite elements, are shown in 6c. of 20 × 20 unit-cells and an excitation of f = 4600 Hz at the center.

3211

It is evident, that both methods provide approximately the same results in the considered frequency range.
However, while the SEM provides the exact results the FEM only delivers approximate solutions. The accuracyof the FEM strongly depends on the number of nodes per wavelength, this leads to an increased computing
time when higher frequency ranges have to be considered.

5

CONCLUSIONSThe properties of a periodic structure can be manipulated by a pre-deformation of the structural elements,
in which the band-gaps can be adjusted by both the amplitude and the shape of the pre-deformation. The
calculation of the corresponding dispersion curves can be done by the Spectral Element Method, which has
the beneﬁt of providing exact results. This is particularly advantageous if higher frequency ranges have to beconsidered.

REFERENCES
1. Xie, L.; et al. An improved fast plane wave expansion method for topology optimization of phononic
crystals. International Journal of Mechanical Sciences, Volume 120, 2016, p. 171-181.
2. Diaz, AR.; Haddow, AG.; Ma, L. Design of band-gap grid structures, Structural and Multidisciplinary Optimization, Vol 29, 2005, p. 418-431.
3. Wang, Y-F.; Wang, Y-S.; Zhang, C. Bandgaps and directional properties of two-dimensional square beamlike zigzag lattices, AIP Advances, Vol 4 (12), 2014, 124403.
4. Yang, C-L.; Zhao, S-D.; Wang, Y-S. Experimental evidence of large complete bandgaps in zig-zag lattice
structures, Ultrasonics, Vol 74, 2017, p. 99-105.
5. Khelif, A.; Adibi, A. Phononic Crystals: Fundamentals and Applications, Springer, New York (USA), 1.
Edition, 2016.
6. Doyle, JF. Wave Propagation in Structures (Chap. 5), Springer, New York (USA), 1. Edition, 1989.
7. Patera, AT. A spectral element method for ﬂuid dynamics - Laminar ﬂow in a channel expansion, Journalof Computational Physics, Vol 54, 1984, p. 468-488.
8. Lee, U. Spectral Element Method in Structural Dynamics, John Wiley & Sons (Asia) Pte Ltd, Singapore(Republic of Singapore), 1. Edition, 2009.
9. Brillouin, L. Wave Propagation in Periodic Structures, McGraw-Hill Book Company, London (United Kingdom), 1. Edition, 1946.
10. Srikantha Phani, A.; Woodhouse J.; Fleck, NA. Wave propagation in two-dimensional periodic lattices,
Acoustical Society of America, Vol 119 (4), 2006, p. 1995-2005.
11. Mellmann, M.; Zhang C. Tuning of vibration and wave propagation characteristics in pre-deformed periodic
lattice frame structures, Proceedings in Applied Mathematics & Mechanics 2019, Vol 19, 2019, Accepted
but not published yet.

3212