arXiv:1409.4812v2 [] 13 Jan 2015

Evaluation of the Spectral Finite Element
Method With the Theory of Phononic
Crystals.
Nicolás Guarı́n-Zapata and Juan Gomez
Departamento de Ingenierı́a Civil
Universidad EAFIT
Medellı́n,
Colombia
jgomezc1@eafit.edu.co
January 15, 2015

Abstract
We evaluated the performance of the classical and spectral finite element
method in the simulation of elastodynamic problems. We used as a quality measure their ability to capture the actual dispersive behavior of the material. Four
different materials are studied: a homogeneous non-dispersive material, a bilayer
material, and composite materials consisting of an aluminum matrix and brass
inclusions or voids. To obtain the dispersion properties, spatial periodicity is
assumed so the analysis is conducted using Floquet-Bloch principles. The effects in the dispersion properties of the lumping process for the mass matrices
resulting from the classical finite element method are also investigated, since that
is a common practice when the problem is solved with explicit time marching
schemes. At high frequencies the predictions with the spectral technique exactly
match the analytical dispersion curves, while the classical method does not. This
occurs even at the same computational demands. At low frequencies however,
the results from both the classical (consistent or mass-lumped) and spectral finite
element coincide with the analytically determined curves.

Keywords: spectral finite elements, numerical dispersion, phononic crystals.
1

1

Introduction

During the recent years, and with particular interest steaming from the earthquake
engineering community, the spectral finite element method (SFEM) has emerged as a
powerful computational tool for numerical simulation of large scale wave propagation
problems (Seriani et al., 1995; Faccioli et al., 1997; Komatitsch et al., 2004; Magnoni
et al., 2012; Cupillard et al., 2012; Shani-Kadmiel et al., 2012; Kudela et al., 2007; Zhu
et al., 2011; Davies et al., 2009; Luo et al., 2009; Luo & Liu, 2009a,b; Balyan et al.,
2012). The origins of the currently used SFEM, can be traced back to the spectral
techniques introduced by Orszag (1980), Patera (1984) and Maday & Patera (1989),
initially proposed for fluid dynamics problems and Gazdag (1981), Kosloff & Baysal
(1982) and Kosloff et al. (1990) in elastodynamics. In seismic wave propagation, the
spectral element methods were first introduced by Priolo et al. (1994) and Faccioli et al.
(1996). Later Komatitsch & Vilotte (1998) replaced the polynomials used by Priolo
et al. (1994) by Legendre polynomials and incorporated a Gauss-Lobatto-Legendre
quadrature, leading to the currently known version of the method in wave propagation
problems. In these methods π-nodal points per wavelength are necessary to resolve
the wave (compared to the 8-10 in the classical FEM) (Ainsworth & Wajid, 2009);
they are also less sensitive to numerical anisotropy and element distortion; exhibit a
smaller conditioning number and lead by construction, to diagonal mass matrices. The
fundamental idea behind the SFEM is the use of higher order Lagrange interpolation at
non-equidistant nodal points, reducing the so-called Runge-phenomenon and producing
exponential convergence rates. If the nodal points, are also the quadrature points
corresponding to a Gauss-Lobatto-Legendre integration rule, the resulting scheme yields
diagonal mass matrices, whereby the global equilibrium equations from explicit time
marching schemes become uncoupled. In contrast, one of the drawbacks of the spectral
technique, is the fact that special meshing methodologies must be developed, since in the
SFEM algorithm the nodal points must be placed at certain (non-standard) positions
within the element. It is thus of interest to identify the range of frequencies for which
the spectral technique results truly advantageous with respect to classical displacement
based formulations. Although there are several works dealing with the performance
of the spectral approach (Marfurt, 1984; Dauksher & Emery, 2000; De Basabe & Sen,
2007; Ainsworth, 2005; Seriani & Oliveira, 2008; Ainsworth & Wajid, 2009; Mazzieri &
Rapetti, 2012), these have been mainly restrained to problems with particular boundary
conditions or else, have been conducted using strongly simplified assumptions.
In this article we evaluate the capabilities of the SFEM in the simulation of mechanical wave propagation problems. We compare its performance with the classical finite
element algorithm using as a quality index, numerically obtained dispersion curves of
the propagation material, computed with the aid of the theory of phononic crystals.
This is an objective performance metric, since it qualifies the ability of the method to
describe the material and not a particular wave propagation problem, thus the effects
2

of geometry and boundary conditions appearing in the solution of a finite model are not
taken into account in this model. For purposes of the evaluation, we implemented classical and spectral finite element algorithms and found numerical dispersion curves for
four different materials, namely a homogeneous material cell, a bilayer material cell, a
porous material cell and a composite material cell. We contrasted the numerical curves
obtained with both methods, with analytic solutions and identified non-dimensional frequency ranges, where the classical and spectral finite element methods deviated from
each other and from the analytic solution. We also evaluated the error, introduced
in the dispersion properties of the material by the artificial diagonalization process
of the finite element mass matrix. This practice is commonly encountered in explicit
time marching schemes implemented with the classical algorithm. From our study we
found that the spectral finite element method yields reliable results at much higher
frequencies than the classical technique even at the same computational costs. At low
frequencies however both methods give the same precision which clearly favours the
classical treatment.
In the first part of the article we briefly describe the implementation of the spectral
finite element method. We then discuss general aspects regarding the numerical formulation of the eigenvalue problem corresponding to the imposition of the Floquet-Bloch
boundary conditions for a spatially periodic material. In the final part of the article we
present results, in terms of dispersion curves for different material cells all of which are
solved numerically and analytically.

2

The Spectral Element Method

The Spectral Finite Element Method (SFEM) is a finite element technique based on high
order Lagrange interpolation (Zienkiewicz & Taylor, 2005) but with a non-equidistant
nodal distribution. The above difference in the spatial sampling, with respect to classical algorithms, seeks to minimize the Runge phenomenon (Pozrikidis, 2005). This
corresponds to the spurious oscillations appearing at the extremes of an interpolation
interval when–in an attempt to perform p-refinement–the number of sampling points is
increased while keeping an equidistant nodal distribution. The problem is illustrated in
Figure 1, where a Runge function is interpolated with 3 different Lagrange polynomials
of order 10, but with a different distribution of the sampling points. We used equidistant nodal points (as in the classical FEM algorithm) and also Lobatto and Chebyshev
nodes (as in the SFEM algorithm).
The Runge phenomenon corresponds to the strong oscillations observed towards
the ends of the interval associated with the equidistant nodal distribution. In contrast,
when the interpolation interval is sampled at the Lobatto and Chebyshev nodes, the
3

approximation describes closely the Runge function throughout the interval. Since
the method uses high-order interpolation for each element, it in fact corresponds to
a p-refinement version of the classical technique, but with a different spatial sampling
inside the elements. The term spectral refers to the rate of convergence since the method
converges faster than any power of 1/p, where p is the order of the polynomial expansion
(Pozrikidis, 2005).
1.5

Runge Function
Equal spaced interpolation
Lobatto Interpolation
Chebyshev Interpolation

y

1

0.5

0

-0.5

-1

-0.5

0

0.5

1

x

Figure 1. Lagrange interpolating polynomials of order 10 for the Runge function
1
.
1+25x2

There are two nodal distributions commonly used in spectral finite element methods: Chebyshev nodes and Lobatto nodes. The Chebyshev nodes are the roots of the
Chebyshev polynomials of the first kind. The resulting interpolation polynomial not
only minimizes the Runge phenomenon but it also provides the best approximating
polynomial under the maximum norm (Burden & Faires, 2011; Kreyzsig, 1989). The
Chebyshev polynomials of the first kind Tn (x) are defined as the unique polynomial
satisfying
Tn (x) = cos [n arccos x]

for each n ≥ 0, and x ∈ [−1, 1] ,

and with roots corresponding to

xk = cos

2k − 1
π
2n


.

Similarly, the Lobatto nodes are the roots of the Lobatto polynomials, which are
the first derivative of the Legendre polynomials
Loi−1 (x) = Pi0 (x) x ∈ [−1, 1] ,
4

with

1 d(i) (t2 − 1)i
Pi (x) = i
,
2 i!
dt(i)
being the ith Legendre Polynomial (Abramowitz & Stegun, 1965). The Lobatto nodes
are also the Fekete points for the line, the square and the cube. The Fekete points
in an interpolation scheme correspond to the nodal distribution that maximizes the
determinant of the Vandermonde matrix (Bos et al., 2001; Pozrikidis, 2005). If the
nodes (corresponding to the sampling points in the interpolation scheme) are also used
as quadrature points in the computation of volume and surface integrals arising in the
FEM, diagonal mass matrices are obtained due to the discrete orthogonality condition
(Pozrikidis, 2005). This last feature of the SFEM makes it very attractive in explicit
time marching schemes, since this condition uncouples the equilibrium equations. In
the classical displacement-based FEM this lumping process is artificially enforced, while
in the SFEM it is a natural result. The nodal distributions corresponding to the 6 × 6
classical and Lobatto-spectral elements are shown in Figure 2.

(a) Classical 6 × 6
element.

(b)
Lobattospectral 6 × 6
element.

Figure 2. Comparison of a classical and a Lobatto-spectral element.

In this work we use the SFEM to solve the reduced elastodynamic frequency domain
equation. The formulation follows from the minimization of the total potential energy
functional Π = Π(ui ) (Reddy, 1991) leading to the principle of virtual displacements
at the frequency ω given by
Z
Z
Z
Z
∗
∗
2
∗
σrs (x)δrs (x)dΩ−ω
ρ(x)ur (x)δur (x)dΩ− tr (x)δur (x)dΓ− fr∗ (x)δur (x)dΩ = 0
Ω

Ω

Γ

Ω

(1)
and where Ω = solution domain, Γ = boundary, σrs = stress tensor, tr = tractions
vector, rs = strain tensor, ur = displacement vector, fr = body force vector and x =
field point vector. The superscripts ∗ appearing in (1) refer to complex conjugate variables. This is a requirement in the weak form, that produces a consistent inner product
between complex variables and leads to self-adjoint operators and therefore Hermitian
finite element matrices: this last condition guarantees the existence of real eigenvalues.
5

The discrete finite element equations follow after writing in (1) the displacements (and
related variables) in terms of basis functions or its derivatives as
ui (x) = NiQ (x)UQ
where ui (x) is the displacement vector evaluated at point x, NiQ (x) is the shape function
for the Qth node and UQ is the nodal displacement vector at the Qth node. In our
notation we retain the index i in the shape function just to indicate that a vectorial
variable is being interpolated. The final finite element equations are of the familiar
form
[K − ω 2 M ] {U} − {F} = {0} .
(2)
where [K] and [M ] are complex valued Hermitian stiffness and mass matrices and F
is the nodal forces vector comprising boundary tractions and body forces. In 2D the
Shape functions for the Q − th node can be computed as the product of independent
Lagrange polynomials like
NiQ (x, y) = `Qx (x)`Qy (y) .

(3)

Similarly, the derivatives of the shape functions are
∂NiQ (x, y)
= `0Qx (x)`Qy (y),
∂x

∂Ni (x, y)
= `Qx (x)`0Qy (y) ,
∂y

(4)

with
`Q (x) =

M
+1
Y
P =1
Q6=P

0Q

` (x) =

M
+1
X
P =1
Q6=P

x − xP
,
xQ − xP

(5)

M
+1
Y
1
x − xR
,
x − xP R=1 xQ − xR

(6)

R6=P

where M is the order of the Lagrange polynomial.
If the integrals in (1) are numerically integrated using the Gauss-Lobatto nodes,
which at the same time are the element nodal points, a diagonal mass matrix is obtained.
Although this is advantageous in the case of an explicit time marching scheme, only
polynomials of order 2N − 3 (Abramowitz & Stegun, 1965; Pozrikidis, 2005) or smaller
are integrated exactly. This accuracy should be contrasted with the one in the standard
Gauss-Legendre quadrature, that integrates exactly polynomials of order 2N − 1 or
smaller (Abramowitz & Stegun, 1965; Bathe, 1995). This reduction in the integration
accuracy is due to the inclusion of nodes in the extremes of the interval. Since the
order of exact integration decreases in the Gauss-Lobatto case, that selection is useful
for interpolations with order higher than N = 3, because the nodal distribution in both
6

methods is the same for N = 1, 2. For the Gauss-Lobatto quadrature the nodes are the
zeros of the completed Lobatto polynomials of order N
(1 − x2 )L0N (x) = 0
with weights
2
N (N + 1)
2
1
wi =
for i = 2, 3, · · · , N .
2
N (N + 1) LN (xi )

w1 = wN +1 =

3

Phononic Crystals

Phononic crystals may be synthetic or natural materials, formed by the spatially periodic repetition of a microstructural feature. Such materials consist of a number of
identical structural components joined together in an identical fashion into an elementary cell that repeats periodically throughout space, forming the complete material and
described by a lattice construction (Mead, 1996). When subjected to waves propagating
at high enough frequencies, the microstructures behave as micro-scatterers, translating into dispersive propagation behaviour at the macroscopic level and sometimes even
leading to states of imaginary group velocity, i.e., values of frequency for which waves do
not propagate denoted as bandgaps. At even higher frequencies, additional propagation
modes may be triggered once a certain cut-off value of the frequency is reached. The
relation between the frequencies at which propagation is possible and the material periodic microstructure–described in terms of the wave vector–can be elegantly described
in terms of the dispersion curve.
In the last few years a strong effort has been devoted to the study of propagation of waves in periodic materials. In particular, methods to determine–via numerical
techniques–the dispersion curves for a specific material, have emerged from earlier work
in the field of solid-state physics, Brillouin (1953). The dispersion curve can be obtained
through the analysis of a single cell, after having identified the microstructure that repeats periodically. Within that context, the unit cell is described by lattice base vectors
indicating the directions that must be followed, in order to cover the complete space
with the application of successive translation operations. As described in Srikantha
et al. (2006), the joints of any lattice structure can be envisioned as a collection of
points (lattice points). These, at the same time are associated with a set of basis vectors ai . The entire direct lattice is then obtained after tessellating the elemental cell
along the basis vectors ai . In this way the location of any point in the periodic material
can be identified with ni integers, defining the number of translations of the elemental
cell along the ai directions of the basis vectors.
7

With the aid of the fundamental tool provided by Bloch’s theorem (Brillouin, 1953)
and extracted from the theory of solid-state physics, the dispersive properties of the
material can be obtained via finite element analysis of the elementary cell. That theorem
states that the change in wave amplitude occurring from cell to cell does not depend
on the wave location within the periodic system (Ruzzene et al., 2003). The final
statement corresponding to Bloch’s theorem, is written in terms of phase shifts between
displacements and tractions on opposite surfaces of the elemental cell like
ui (x + a) = ui (x)eik·a ,
ti (x + a) = −ti (x)eik·a

(7)

and where k is the wave vector and a = a1 n1 + a2 n2 + a3 n3 is the lattice translation
vector, formed after scaling the basis vectors by the translation integers ni . The factor
eik·a corresponds to the phase shift between opposite sides of the cell. The real and
imaginary components of the wave vector differ by this factor, while the magnitude
remains the same.

3.1

Bloch-periodicity in discrete methods

The formulation of the finite element equations for the determination of the dispersion
properties of an elemental cell requires the assembly of the discrete equations for the
complete cell, the subsequent reduction of the equations after imposition of the Blochperiodic boundary conditions stated in (7), and the solution of a series of generalized
eigenvalue problems to find the dispersion curve. If the eigenvalue problem or numerical
dispersion relation is written like
D̂(k, ω)U = 0

(8)

where k is the wave vector which is progressively assigned successive values, in such a
way that the first Brillouin zone is fully covered. Each application of the wave vector
and the solution of the related eigenvalue problem yields a duple (k, ω) representing a
plane wave propagating at frequency ω.
The basic steps are further elaborated with reference to Figure 3, describing a unit
cell and where 2da and 2db are respectively the cell width and cell height. The subindices labelled as 9 (or i) are related to the internal degrees of freedom, while those
labelled 1 − 4 ( or alternatively l, r, b, t) are used for degrees of freedom over the boundary and 5 − 8 (or lb, rb, lt, rt) correspond to those over the corners of the unit cell
(Srikantha et al., 2006). In the initial step we assembly the complete discrete finite element equations for the unit cell using the following block vectors of nodal displacements

8

and forces
U = [ul ur ub ut ulb urb ult urt ui ]T ≡ [u1 u2 u3 u4 u5 u6 u7 u8 u9 ]T ,
F = [fl fr fb ft flb frb flt frt fi ]T ≡ [f1 f2 f3 f4 f5 f6 f7 f8 f9 ]T .

ult u7

ut u4

urt u8

ult u7

ut u4

urt u8

ul u1

ui u9

ur u2

ul u1

ui u9

ur u2

urb u6

ulb u5

ulb u5

ub u3

(a) Original cell.

ult
ul

ulb

ub u3

urb u6

(b) Reduced cell.

u

21

22

23 t 24

25

16

17

18

19

20

14

15
10

ui

11

12

13

6

7

8

9

1

2

3

ub 4

5

urt
ur

urb

(c) The black circles now correspond to single degrees of freedom.
Notice that the nodal numbering does not match the original numbering for groups of degrees of freedom.
Figure 3. Relevant sets of DOF for an schematic unitary cell discretized with several
finite elements. Each circle represents a family of degrees of freedom for a typical
mesh and not a single degree of freedom. In part a) we show all the degrees of
freedom for the unitary cell before imposing the relevant Bloch-periodic boundary
conditions. Part b) shows the reduced cell where the white circles enclosed by the
dark square represent the reference nodes containing the information from the image
nodes which will be eventually deleted from the system. Part c) shows an example
of node grouping for a mesh of 4 × 4 elements.

In the FEM equations the Bloch conditions take the form of constraint relationships
among the degrees of freedom in opposite boundaries of the cell. This allows reduction
of the system of equations via elimination of the subset of image nodes and leaving only

9

the subset of reference nodes. These constraints are written like
u2 = eiψx u1
u4 = eiψy u3
u6 = eiψx u5
u7 = eiψy u5
u8 = ei(ψx +ψy ) u5 ,
where ψx = 2kx da and ψy = 2ky db define phase shifts in x and y directions respectively
and (kx , ky ) = k corresponds to the wave vector. Similarly, using the corresponding
relationships between the nodal forces we have
f2 + eiψx f1 = 0
f4 + eiψy f3 = 0
f8 + eiψx f6 + eiψy f7 + ei(ψx +ψy ) f5 = 0.
If we assume that the problem does not involve internal forces, i.e. f9 = 0, we get
the reduced generalized eigenvalue problem to be solved for each frequency and related
propagation mode and given by
[KR − ω 2 MR ]{UR } = {0} .

(9)

In a finite element implementation the Bloch-periodic boundary conditions can be considered either including directly the phase shifts in the basis functions at the outset
or performing row and column operations in the assembled matrices, Sukumar & Pask
(2009). The Hermitian matrices KR and MR in the reduced system are function of the
wave vector and angular frequency as indicated in (8).

4

Results

In order to test the performance of the spectral finite element method, we implemented
an in-house code to assemble and solve the generalized eigenvalue problem stated in
(9). We considered a homogeneous (non-dispersive material), a bilayer material and
a composite material made out of an aluminium matrix with a square inclusion. In
particular we studied the case of a pore and a brass inclusion. From the above set, the
homogeneous and bilayer material cells had closed form solutions, so we were able to
compare our numerical results with the analytical dispersion curves.
In all the considered cases we solved the problem using both the classical finite
element technique and the spectral finite element method. On the other hand, and
10

since one of the reported advantages of the spectral technique, with respect to the
classical approach, is the fact that in this method a diagonal mass matrix is obtained,
we also tested the effect of using artificially lumped mass matrices in the classical FEM.
With this goal in mind, we also solved the homogeneous cell and the bilayer material
cell, using the classical FEM with a consistent mass matrix and a lumped mass matrix.
The method for making the matrices diagonals is the Diagonal scaling (see Zienkiewicz
et al. (2005)), where each of the terms over the diagonal are scaled by a factor under
the constraint of keep the total mass of the element constant. In our case, the factor is
the total mass of the element (Mtot ) divided by the trace of the matrix (Tr(M )), i.e.
(lumped)

Mii

=

Mtot
Mii
Tr(M )

(no summation on i) .

The results are presented next in terms of dimensionless frequency vs dimensionless
wave number.

4.1

Homogeneous material

In an ideal homogeneous material a plane propagates non-dispersively according to the
linear relation
ω = ckkk.
(10)
If the wave is travelling parallel to the y-axis (i.e., vertical incidence) and we consider
a square cell of side 2d the above relation particularizes to
r
nπ 2  mπ 2
k+
+
,
(11)
ω i = ci
d
d
where ci is the wave propagation velocity for the P wave (i = 1) or SV wave (i = 2). In
this relation n and m are integers indicating the relative position of any given cell with
respect the reference cell. When n is different from zero, the dispersion branches become
hyperbolas due to the wave vector whose components are determined as modulus 2π and
could be present in any Brillouin zone. These hyperbolas result from the intersection
of a plane parallel to the x−axis and a cone, Langlet (1993); Langlet et al. (1995);
Guarı́n-Zapata (2012).
Figure 4 presents the results for a square cell discretized with 4-th and 7-th order
complete Lagrange polynomials with equidistant and Lobatto sampling nodes. The
analytical and numerical dispersion curves are shown in dotted and continuous lines
respectively. In all the considered cases the differences between the numerical and analytical results are due to numerical dispersion. The results from the classical 5×5
11

element are accurate only at normalized frequencies below 3.0, while the analogous
analysis with the spectral method predicts even modes of higher multiplicity as seen
from the point near (1,4). Moreover, the classical 8×8 element only captures the first
three modes, while this same result was already reached with the 5×5 spectral element.
This of course is a result of the improved convergence properties of the spectral technique. Similarly the 8×8 spectral element captures the solution exactly up to the sixth
mode.

6

6

5

5

4

4

3

3

2

2

1

1

0
0.0

0.2

0.4

0.6

0.8

1.0

(a) Classical 5 × 5 element.

0
0.0

6

5

5

4

4

3

3

2

2

1

1

0.2

0.4

0.6

0.4

0.6

0.8

1.0

(b) Spectral 5 × 5 element.

6

0
0.0

0.2

0.8

1.0

(c) Classical 8 × 8 element.

0
0.0

0.2

0.4

0.6

0.8

(d) Spectral 8 × 8 element.

Figure 4. Dispersion curves from a single homogeneous material cell.

12

1.0

4.2

Bilayer material

Figure 5 describes a composite material formed by two layers of different properties.
When the material cell is submitted to an incident wave in the direction perpendicular to
the layers the dispersion curve can be found in closed-form. The problem was studied by
Lord Rayleigh in the case of electromagnetic waves which, in this case, is mathematically
equivalent to the mechanical problem (Rayleigh, 1888). Of particular interest in this

Figure 5. Cell for a bilayer material.

material is the presence of band gaps in addition of its dispersive behaviour. The
analytic dispersion relation as obtained by Rayleigh (1888) is given by
 
 
 
 
ωd
(ρ1 c1 )2 + (ρ2 c2 )2
ωd
ωd
ωd
cos
−
sin
sin
,
cos(2dk) = cos
c1
c2
2ρ1 ρ2 c1 c2
c1
c2
where ci the longitudinal/transversal wave speed of the i layer and ρi the density of
the i layer. In the current analysis the material properties are EA = 7.31 × 1010 Pa,
νA = 0.325, ρA = 2770 kg/m3 and EB = 9.2 × 1010 Pa, νB = 0.33, ρB = 8270
kg/m3 for aluminium and brass respectively. Figure 6 shows the comparison between
the results obtained with classical FEM (with consistent and lumped mass matrix)
SFEM and the analytical solution. Notice that the same number of modes (5 modes),
were found with the classical 8 × 8 finite element and with the lower order 4 × 4
spectral element. Furthermore, the spectral element, as shown in Figure 7, reproduces
exactly up to 11 modes, once again showing the improved convergency properties of
the spectral method. It is evident that in the considered cases this approximation does
not introduce considerable error in the dispersive properties as compared with the full
consistent matrix.

13

6

6

5

5

4

4

3

3

2

2

1

1

0
0.0

0.2

0.4

0.6

0.8

1.0

(a) Classical 4 × 4 element.

0
0.0

6

5

5

4

4

3

3

2

2

1

1
0.2

0.4

0.6

0.8

0.4

0.6

0.8

1.0

(b) Classical 8 × 8 element.

6

0
0.0

0.2

1.0

0
0.0

0.2

0.4

0.6

0.8

1.0

(c) Classical 4 × 4 lumped element. (d) Classical 8 × 8 lumped element.
6

6

5

5

4

4

3

3

2

2

1

1

0
0.0

0.2

0.4

0.6

0.8

1.0

(e) Spectral 4 × 4 element.

0
0.0

0.2

0.4

0.6

0.8

1.0

(f ) Spectral 8 × 8 element.

Figure 6. Dispersion curves for a bilayer material. Each layer is meshed with a
single element with homogeneous material, being the materials aluminium and brass.

14

12
10
8
6
4
2
0
0.0

0.2

0.4

0.6

0.8

1.0

Figure 7. Results for 12 branches of the dispersion relation for the bilayer material
computed with a single 8 × 8 SFEM element for layer. The materials are aluminium
and brass.

4.3

Matrix with pore/elastic inclusion

As a final case we considered the cell of an aluminium matrix with a square pore and
a square brass inclusion. Figure 8 shows the dispersion curves with the 8 × 8 spectral
element used to capture the 11 modes in the previous case. This result is also reported
in Langlet (1993). There is a bandgap with a cut-off frequency close to 3.0 in the cell
with the pore. The bandgap however, completely closes in the case of the inclusion
since it presents a smaller impedance contrast.

5

5

4

4

3

3

2

2

1

1

0
0.0

0.2

0.4

0.6

0.8

0
0.0

1.0

(a) Pore.

0.2

0.4

0.6

0.8

1.0

(b) Inclusion.

Figure 8. Dispersion relations for a square brass inclusion and pore in a square cell
of aluminium. Both were computed with 8 × 8 SFEM elements.

15

5

Conclusions

We have implemented a higher order classical finite element method and spectral finite element method to find the dispersion relations for in-plane P and SV waves,
propagating in periodic materials, i.e, phononic crystals. The implemented solution
method is based on the Floquet-Bloch theorem from solid state physics, which allows
the determination of the dispersion curve from the analysis of a single cell. The numerical dispersion curves of the classical and spectral approaches were used to qualify
the performance of each technique in the simulation of mechanical wave propagation
problems. Since the dispersion curve simultaneously describes properties of the propagation medium and the external excitation, it is an engineering and objective measure
to assess the performance of the numerical method and particularly of its ability to
correctly capture the dispersive properties of different material microstructures.
The considered materials correspond to a homogeneous non-dispersive material, a
dispersive bilayer material and a composite material exhibiting bandgaps. The composite material is formed by a matrix with a square pore/inclusion. In general, it was
found that the spectral finite element method yields reliable results at much higher
frequencies than the classical technique, even at the same computational costs. At low
frequencies however, the performance of both methods is equivalent and the analysis
could rely on the classical algorithm. One of the main advantages of the SFEM is
the fact that having used Gauss-Lobatto-Legendre nodes the resulting mass matrix is
diagonal by formulation. In the case of an analysis using Bloch-periodic boundary conditions, a diagonal mass matrix allows reducing the generalized eigenvalue problem into
a standard eigenvalue problem. This can be accomplished via Cholesky decomposition.
As a secondary goal, we also evaluated in this work the effects of enforcing diagonality of the mass matrix resulting from a classical algorithm. The dispersion curves
obtained with the classical method using a consistent and a lumped mass matrix were
the same in practical terms which allows to use such approximation with confidence.

6

Acknowledgements

This work was supported by Colciencias and Universidad EAFIT under the Young
Researcher and Innovators Program “Virginia Gutiérrez de Pineda”. The authors would
like to thank Professor Anne Christine Hladky-Hennion from Institut d’électronique de
microélectronique et de nanotechnologie -Université Lille Nord, France for answering so
kindly our conceptual queries and also providing us a copy of Langlet (1993).

16

References
Abramowitz, M. & Stegun, I., 1965. Handbook of mathematical functions: with formulas, graphs, and mathematical tables, vol. 55, Dover publications.
Ainsworth, M., 2005. Discrete dispersion relation for hp-version finite element approximation at high wave number, SIAM Journal on Numerical Analysis, 42(2), 553–575.
Ainsworth, M. & Wajid, H., 2009. Dispersive and dissipative behavior of the spectral
element method, SIAM Journal on Numerical Analysis, 47(5), 3910–3937.
Balyan, L. K., Dutt, P. K., & Rathore, R., 2012. Least squares h-p spectral element method for elliptic eigenvalue problems, Applied Mathematics and Computation, 218(19), 9596 – 9613.
Bathe, K.-J., 1995. Finite Element Procedures, Prentice Hall, 2nd edn.
Bos, L., Taylor, M., & Wingate, B., 2001. Tensor product Gauss-Lobatto points are
Fekete points for the cube, Mathematics of Computation, 70(236), 1543–1548.
Brillouin, L., 1953. Wave propagation in periodic structures: Electric filters and crystal
lattices, Dover Publications, 1st edn.
Burden, R. & Faires, J., 2011. Numerical Analysis, Brooks/Cole, 9th edn.
Cupillard, P., Delavaud, E., Burgos, G., Festa, G., Vilotte, J.-P., Capdeville, Y., & Montagner, J.-P., 2012. RegSEM: a versatile code based on the spectral element method
to compute seismic wave propagation at the regional scale, Geophysical Journal International , 188(3), 1203–1220.
Dauksher, W. & Emery, A., 2000. The solution of elastostatic and elastodynamic problems with chebyshev spectral finite elements, Computer methods in applied mechanics
and engineering, 188(1), 217–233.
Davies, R., Morgan, K., & Hassan, O., 2009. A high order hybrid finite element method
applied to the solution of electromagnetic wave scattering problems in the time domain, Computational Mechanics, 44, 321–331.
De Basabe, J. D. & Sen, M. K., 2007. Grid dispersion and stability criteria of some
common finite-element methods for acoustic and elastic wave equations, Geophysics,
72(6), T81–T95.
Faccioli, E., Maggio, F., Quarteroni, A., & Tagliani, A., 1996. Spectral-domain decomposition methods for the solution of acoustic and elastic wave equations, Geophysics,
61(4), 1160–1174.

17

Faccioli, E., Maggio, F., Paolucci, R., & Quarteroni, A., 1997. 2d and 3d elastic
wave propagation by a pseudo-spectral domain decomposition method, Journal of
seismology, 1(3), 237–251.
Gazdag, J., 1981. Modeling of the acoustic wave equation with transform methods,
Geophysics, 46(6), 854–859.
Guarı́n-Zapata, N., 2012. Simulación Numérica de Problemas de Propagación de Ondas:
Dominios Infinitos y Semi-infinitos, Master’s thesis, Universidad EAFIT.
Komatitsch, D. & Vilotte, J., 1998. The spectral element method: An efficient tool
to simulate the seismic response of 2d and 3d geological structures, Bulletin of the
Seismological Society of America, 88(2), 368–392.
Komatitsch, D., Liu, Q., Tromp, J., Süss, P., Stidham, C., & Shaw, J. H., 2004.
Simulations of Ground Motion in the Los Angeles Basin Based upon the SpectralElement Method, Bulletin of the Seismological Society of America,, 94(1), 187–206.
Kosloff, D. & Baysal, E., 1982. Forward modeling by a fourier method, Geophysics,
47(10), 1402–1412.
Kosloff, D., Kessler, D., Tessmer, E., Behle, A., Strahilevitz, R., et al., 1990. Solution
of the equations of dynamic elasticity by a chebychev spectral method, Geophysics,
55(6), 734–747.
Kreyzsig, E., 1989. Introductory Functional Analysis with Applications, Wiley, 1st edn.
Kudela, P., Zak, A., Krawczuk, M., & Ostachowicz, W., 2007. Modelling of wave
propagation in composite plates using the time domain spectral element method,
Journal of Sound and Vibration, 302(4-5), 728 – 745.
Langlet, P., 1993. Analyse de la Propagation des Ondes Acoustiques dans les Materiaux
Periodiques a l’aide de la Methode des Elements Finis, Ph.D. thesis, L’Universite de
Valenciennes et du Hainaut-Cambresis.
Langlet, P., Hladky-Hennion, A.-C., & Decarpigny, J.-N., 1995. Analysis of the propagation of plane acoustic waves in passive periodic materials using the finite element
method, Journal of the Acoustical Society of America, 98(5), 2792–2800.
Luo, M. & Liu, Q. H., 2009a. Accurate determination of band structures of twodimensional dispersive anisotropic photonic crystals by the spectral element method,
Journal of The Optical Society of America, 26(7), 1598–1605.
Luo, M. & Liu, Q. H., 2009b. Spectral element method for band structures of threedimensional anisotropic photonic crystals, Physical Review E , 80, 056702.
Luo, M., Liu, Q. H., & Li, Z., 2009. Spectral element method for band structures of
two-dimensional anisotropic photonic crystals, Physical Review E , 79, 026705.
18

Maday, Y. & Patera, A. T., 1989. Spectral element methods for the incompressible Navier-Stokes equations, in State-of-the-art surveys on computational mechanics
(A90-47176 21-64). New York, American Society of Mechanical Engineers, 1989, p.
71-143. Research supported by DARPA., pp. 71–143.
Magnoni, F., Casarotti, E., Michelini, A., Piersanti, A., Komatitsch, D., Peter, D., &
Tromp, J., 2012. Spectral-element and adjoint 3d tomography for the lithosphere of
central italy, in EGU General Assembly Conference Abstracts, vol. 14, p. 7506.
Marfurt, K. J., 1984. Accuracy of finite-difference and finite-element modeling of the
scalar and elastic wave equations, Geophysics, 49(5), 533–549.
Mazzieri, I. & Rapetti, F., 2012. Dispersion analysis of triangle-based spectral element
methods for elastic wave propagation, Numerical Algorithms, 60(4), 631–650.
Mead, D., 1996. Wave propagation in continuous periodic structures: research contributions from southampton, 1964–1995, Journal of Sound and Vibration, 190(3),
495–524.
Orszag, S., 1980. Spectral methods for problems in complex geometries, Journal of
Computational Physics, 37(1), 70–92.
Patera, A. T., 1984. A spectral element method for fluid dynamics: Laminar flow in a
channel expansion, Journal of Computational Physics, 54, 468 – 488, This was the
first paper in Spectral Element Methods.
Pozrikidis, C., 2005. Introduction To Finite And Spectral Element Methods Using Matlab, Chapman & Hall/CRC.
Priolo, E., Carcione, J., & Seriani, G., 1994. Numerical simulation of interface waves
by high-order spectral modeling techniques, The Journal of the Acoustical Society of
America, 95, 681.
Rayleigh, L., 1888. On the remarkable phenomenon of crystalline reflexion described
by Prof. Stokes, Philosophical Magazine Series 5 , 26(160), 256–265.
Reddy, J., 1991. Applied Functional Analysis and Variational Methods in Engineering,
Krieger Publishing, 1st edn.
Ruzzene, M., Scarpa, F., & Soranna, F., 2003. Wave beaming effects in two-dimensional
cellular structures, Smart Materials and Structures, 12, 363–372.
Seriani, G. & Oliveira, S., 2008. Dispersion analysis of spectral element methods for
elastic wave propagation, Wave Motion, 45(6), 729–744.
Seriani, G., Priolo, E., & Pregarz, A., 1995. Modelling waves in anisotropic media by
a spectral element method, in Proceedings of the third international conference on
mathematical and numerical aspects of wave propagation, pp. 289–298, Siam.
19

Shani-Kadmiel, S., Tsesarsky, M., Louie, J., & Gvirtzman, Z., 2012. Simulation of
seismic-wave propagation through geometrically complex basins: The dead sea basin,
Bulletin of the Seismological Society of America, 102(4), 1729–1739.
Srikantha, A., Woodhouse, J., & Flecka, N., 2006. Wave propagation in two-dimensional
periodic lattices, Journal of the Acoustical Society of America, 119, 1995–2005.
Sukumar, N. & Pask, J. E., 2009. Classical and enriched Finite element formulations
for Bloch-periodic boundary conditions, International Journal of Numerical Methods
in Engineering, 77(8), 1121–1138.
Zhu, C., Qin, G., & Zhang, J., 2011. Implicit chebyshev spectral element method for
acoustics wave equations, Finite Elements in Analysis and Design, 47(2), 184 – 194.
Zienkiewicz, O. & Taylor, R., 2005. The Finite Element Method for Solid and Structural
Mechanics, vol. 2, Butterworth-Heinemann, 6th edn.
Zienkiewicz, O., Taylor, R., & Zhu, J., 2005. The Finite Element Method: Its Basis
and Fundamentals, vol. 1, Butterworth-Heinemann, 6th edn.

20

