INTERNATIONAL JOURNAL FOR NUMERICAL METHODS IN ENGINEERING
Int. J. Numer. Meth. Engng 2010; 00:1–26
Published online in Wiley InterScience (www.interscience.wiley.com). DOI: 10.1002/nme
A Discontinuous Galerkin Method with Plane Waves
for Sound Absorbing Materials
G. Gabard
1
, O. Dazel
2
1
ISVR, University of Southampton, Southampton, UK
2
LAUM, UMR CNRS 6613, Universit
´
e du Maine, Le Mans, France
SUMMARY
Poro-elastic materials are commonly used for passive control of noise and vibration, and are key to
reducing noise emissions in many engineering applications, including the aerospace, automotive and energy
industries. More efficient computational models are required to further optimise the use of such materials. In
this paper we present a Discontinuous Galerkin method (DGM) with plane waves for poro-elastic materials
using the Biot theory solved in the frequency domain. This approach offers significant gains in computational
efficiency and is simple to implement (costly numerical quadratures of highly-oscillatory integrals are not
needed). It is shown that the Biot equations can be easily cast as a set of conservation equations suitable for
the formulation of the wave-based DGM. A key contribution is a general formulation of boundary conditions
as well as coupling conditions between different propagation media. This is particularly important when
modelling porous materials as they are generally coupled with other media, such as the surround fluid or
an elastic structure. The validation of the method is described first for a simple wave propagating through
a porous material, and then for the scattering of an acoustic wave by a porous cylinder. The accuracy,
conditioning and computational cost of the method are assessed, and comparison with the standard finite
element method is included. It is found that the benefits of the wave-based DGM are fully realised for the
Biot equations and that the numerical model is able to accurately capture both the oscillations and the rapid
attenuation of the waves in the porous material. Copyright
c
2010 John Wiley & Sons, Ltd.
Received . . .
KEY WORDS: porous material, Biot theory, discontinuous Galerkin method, plane wave
1. INTRODUCTION
The objective of this work is to develop a Discontinuous Galerkin Method (DGM) with plane waves
to predict sound absorption in poro-elastic materials (PEM). Such materials are commonly used
for passive control of noise and vibration. In practical applications, they are often combined as
layers attached or linked to a vibrating structure or to a fluid cavity. When subjected to mechanical
or acoustical excitation, they can dissipate energy through viscous, thermal and structural effects,
making their computational modelling of particular importance for many engineering applications.
The dynamic behaviour of porous materials is classically obtained from homogenized models and
particularly from the Biot–Allard theory [1–4] which is based on a continuous field mechanics
approach. The homogenized porous media is modelled as the combination of two continuous fields
whose inertial and constitutive coefficients are given by phenomenological relations.
Prediction models for porous materials, by analogy with structures, are commonly classified
according to three categories. The first corresponds to low frequencies for which the main
computational technique is the Finite Element Method (FEM). The third category is associated
∗
Correspondence to: Email: olivier.dazel@univ-lemans.fr
Copyright
c
2010 John Wiley & Sons, Ltd.
Prepared using nmeauth.cls [Version: 2010/05/13 v3.00]
2 G. GABARD, O. DAZEL
with high frequencies where the Transfer Matrix Method provides an efficient representation of the
layers of materials mentioned above. Standing between these two extremes, the second category is
commonly referred to as the mid-frequency range and is currently the subject of active research.
The boundaries between these regimes are somewhat arbitrary, as they depend on the properties
and geometry of the porous material and of the structure it is attached to. The mid-frequencies
correspond to situations where the standard FEM struggles with the size of the problem, but the
statistical methods commonly used for high-frequency problems are not yet applicable. Most of the
computational methods for porous materials discussed in the literature are extensions of either the
low- or high-frequency methods. The method proposed in the present paper is designed for the low-
and mid-frequency regimes.
Discontinuous Galerkin Methods (DGM) have been actively developed for various branches
of science, especially for time-domain simulations of conservation equations, as these methods
directly provide high-order, explicit schemes [5]. They generally rely on polynomial interpolations
of the solution within each element, but recently the use of plane waves has been proposed.
This is part of the development of so-called wave-based, or Trefftz, methods where the use of
canonical solutions as basis functions improves significantly the accuracy of the numerical model,
and offers exponential convergence when the number of plane waves is increased. Prominent
examples of this approach include the partition of unity finite element method (PUFEM) [6],
the ultra-weak variational formulation [7] and the discontinuous enrichment method [8]. The use
of the discontinuous Galerkin approach in this context was then proposed using either Lagrange
multipliers [9] or numerical flux methods [10]. More recently a thorough analysis of the properties
of the wave-based DGM for the Helmholtz equation was conducted [11–14]. It was also shown
that the DGM with numerical flux provides a unified framework to describe several wave-based
methods [10, 15, 16], including the ultra-weak variational formulation, and the wave-based least-
square method [17].
In this paper we present a DGM using plane waves for poro-elastic materials modelled with
Biot’s theory in the frequency domain. Relevant prior work include the use of a Trefftz wave-
based approach [18, 19] as well as the ultra-weak formulation for the lossy wave equation [20].
The PUFEM was recently applied to an equivalent fluid model [21] and then to the full Biot
equations [22]. An issue with the PUFEM is the cost of calculating the element matrices associated
with numerical quadratures for highly-oscillatory integrands (it was reported in some cases that
the cost of calculating the element matrices is of the same order as solving the system of linear
equations). The present wave-based DGM does not suffer from this issue.
The Biot equations are introduced in the next section and it is shown that they can easily be cast
as a set of conservation equations which is needed for the present DGM. In section 3 the general
formulation of the wave-based DGM is recalled, and the emphasis is placed on the formulation
of boundary conditions as well as coupling conditions between two different propagation media
(typically between the porous medium and the surround fluid). The validation of the method is
described in section 4, first for a simple wave propagating through the PEM, and then for the
scattering of an acoustic wave by a porous cylinder. The accuracy, conditioning and computational
cost of the method are assessed.
2. GOVERNING EQUATIONS
Throughout this paper, a harmonic time dependence e
+iωt
is assumed with the angular frequency
ω. The numerical method and its applications are presented in two dimensions (x, y). To model the
propagation of waves in the porous material and in the surrounding fluid, we will consider a general
system of linear conservation equations of the form:
iωu + A
x
∂u
∂x
+ A
y
∂u
∂y
= 0 , (1)
where u is the vector of physical field variables whose number depends on the nature of the
propagation medium. Below we consider either a porous material, or a fluid (using 8 or 3 field
Copyright
c
2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010)
Prepared using nmeauth.cls DOI: 10.1002/nme
PLANE-WAVE DGM FOR POROELASTIC MATERIALS 3
variables, respectively). The coefficient matrices A
x
and A
y
are assumed constant. For the porous
material, these matrices are complex valued and vary with frequency ω. The conservation equations
(1) represent the basis for most discontinuous Galerkin methods, and we will also use this as a
starting point to formulate the proposed wave-based DGM in section 3.
2.1. First-order model for poro-elastic media
There are several formulations of the Biot equations available in the literature. In the present work
we will use the Biot equations as formulated in [23] as the simplicity of this formulation greatly
facilitates the algebra when deriving the plane wave basis. The only difference with [23] is that the
velocity is used here instead of the displacement. The equations of motion involve the velocity v
s
of the solid phase of the porous material, as well as the total velocity v
t
of the porous material:
iωeρ
s
v
s
x
+ iωeρ
eq
eγv
t
x
=
∂ˆσ
xx
∂x
+
∂ˆσ
xy
∂y
, (2a)
iωeρ
s
v
s
y
+ iωeρ
eq
eγv
t
y
=
∂ˆσ
xy
∂x
+
∂ˆσ
yy
∂y
, (2b)
iωeγ eρ
eq
v
s
x
+ iωeρ
eq
v
t
x
= −
∂p
f
∂x
, (2c)
iωeγ eρ
eq
v
s
y
+ iωeρ
eq
v
t
y
= −
∂p
f
∂y
. (2d)
The left-hand side of these equations is related to visco-inertial terms. The PEM equivalent densities
˜ρ
s
, ˜ρ
eq
and the coupling factor ˜γ are defined in [23]. The right-hand side corresponds to elastic
effects. The pressure in the fluid phase is denoted p
f
. The tensor
ˆ
σ corresponds to the stresses in
the solid phase of the porous material in the absence of fluid (i.e. in vacuo). These are defined as
follows:
iωp
f
= −
e
K
eq
∂v
t
x
∂x
+
∂v
t
y
∂y
, (3a)
iωˆσ
xx
=
ˆ
A
∂v
s
y
∂y
+ (
ˆ
A + 2N)
∂v
s
x
∂x
, (3b)
iωˆσ
xy
= N
∂v
s
x
∂y
+
∂v
s
y
∂x
, (3c)
iωˆσ
yy
= (
ˆ
A + 2N)
∂v
s
y
∂y
+
ˆ
A
∂v
s
x
∂x
. (3d)
e
K
eq
is the compressibility of the fluid.
ˆ
A and N are the elastic coefficients of the solid phase in vacuo
which are directly obtained from the Lam
´
e coefficients λ and µ, as shown in [23] and recalled in
Appendix A.
From equations (2) and (3) it is clear that can first be cast into a non-conservative form
iωMu + B
x
∂u
∂x
+ B
y
∂u
∂y
= 0 , (4)
with the following definitions:
u =
v
s
x
v
s
y
v
t
x
v
t
y
ˆσ
+
ˆσ
xy
ˆσ
−
p
f
, M =
eρ
s
0 eγeρ
eq
0 0 0 0 0
0 eρ
s
0 eγeρ
eq
0 0 0 0
eγeρ
eq
0 eρ
eq
0 0 0 0 0
0 eγeρ
eq
0 eρ
eq
0 0 0 0
0 0 0 0 (
ˆ
A + N)
−1
0 0 0
0 0 0 0 0 N
−1
0 0
0 0 0 0 0 0 N
−1
0
0 0 0 0 0 0 0
e
K
−1
eq
, (5)
Copyright
c
2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010)
Prepared using nmeauth.cls DOI: 10.1002/nme
4 G. GABARD, O. DAZEL
and
B
x
=
0 0 0 0 −1 0 −1 0
0 0 0 0 0 −1 0 0
0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0
−1 0 0 0 0 0 0 0
0 −1 0 0 0 0 0 0
−1 0 0 0 0 0 0 0
0 0 1 0 0 0 0 0
, B
y
=
0 0 0 0 0 −1 0 0
0 0 0 0 −1 0 1 0
0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1
0 −1 0 0 0 0 0 0
−1 0 0 0 0 0 0 0
0 1 0 0 0 0 0 0
0 0 0 1 0 0 0 0
.
(6)
The reason for introducing the non-conservative form (4) is that deriving closed-form expression
of the plane-wave basis for the solution and the test functions is easier when using (4), since it
provides simple links between the direct and adjoint plane wave bases, as will be shown in section
3.6. This stems from the fact that the matrices B
x
and B
y
are real and symmetric, while the complex-
valued matrix M is complex symmetric (but not Hermitian).
The conservative form (1) is easily recovered from equation (4) using
A
x
= M
−1
B
x
, A
y
= M
−1
B
y
. (7)
For convenience we have introduced σ
+
= (ˆσ
xx
+ ˆσ
yy
)/2 and σ
−
= (ˆσ
xx
− ˆσ
yy
)/2 in the vector
u as it simplifies the expression of the mass matrix M and the calculation of its inverse in (7).
It should be noted that it is also possible to formulate a DGM starting from the non-conservative
form (4), provided that the numerical flux is defined in a consistent way to ensure conservation
of the field variables. This approach has been derived and implemented by the authors, but it is
not described in the present paper as this leads to weak forms equivalent to that obtained from the
well-established DGM based on (1).
One might think that the large number of unknowns introduced in this model implies that the
computational cost of solving for all these variables will be high. This is true for standard finite
element methods where each variable is discretised independently, but it does not apply here. With
the present wave-based method the degrees of freedom are the amplitudes of the plane waves in
each element and their number is completely independent of the number of variables introduced in
the governing equations.
2.2. Acoustic waves in air
To describe the acoustic waves in the fluid around the porous material we use the standard Helmholtz
equation which can be written directly in the conservative form (1) by introducing the acoustic
pressure p
a
and linearised momentum ρ
0
v
a
as field variables:
u =
p
a
ρ
0
v
a
x
ρ
0
v
a
y
, A
x
=
0 c
2
0
0
1 0 0
0 0 0
, A
y
=
0 0 c
2
0
0 0 0
1 0 0
, (8)
where ρ
0
is the mean density and c
0
is the sound speed. This corresponds to the same set of equations
used in [16].
3. WAVE-BASED DGM
We will now present the formulation and discretisation of the wave-based discontinuous Galerkin
method of the conservative equations (1). We follow the same principles as in [10, 16], but the
significant addition presented here is a general approach to incorporate a large class of boundary
conditions (section 3.5), as well as coupling conditions between two different media (section 3.4).
This approach relies heavily on the concept of characteristics which is introduced in section 3.2.
Copyright
c
2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010)
Prepared using nmeauth.cls DOI: 10.1002/nme
PLANE-WAVE DGM FOR POROELASTIC MATERIALS 5
3.1. Variational formulation
We consider a domain Ω which is represented by a set of N
e
elements Ω
e
. We allow for the solution
u to be discontinuous at the interfaces between the elements. The variational formulation associated
with the conservative form (1) is to find a solution u such that
X
e
Z
Ω
e
v
T
e
iωu
e
+ A
x
∂u
e
∂x
+ A
y
∂u
e
∂y
dΩ = 0 , ∀v , (9)
where
T
denotes the Hermitian transpose. u
e
= u|
Ω
e
and v
e
= v|
Ω
e
denote the restrictions of the
solution and the test function to each element Ω
e
.
After integrating by parts on each element and rearranging terms we get:
−
X
e
Z
Ω
e
iωv
e
+ A
T
x
∂v
e
∂x
+ A
T
y
∂v
e
∂y
T
u
e
dΩ +
X
e
Z
∂Ω
e
v
T
e
F
e
u
e
dΓ = 0 , ∀v , (10)
where we have introduced the matrix F
e
= A
x
n
x
+ A
y
n
y
which represents the normal fluxes
across the boundary of the element Ω
e
. The unit normal n = (n
x
, n
y
) on the element boundary
∂Ω
e
points out of the element.
A key aspect of the wave-based DGM is to use test functions v whose restrictions v
e
on each
elements are solutions of the adjoint problem defined on each element:
iωv
e
+ A
T
x
∂v
e
∂x
+ A
T
y
∂v
e
∂y
= 0 , (11)
which is readily identified from equation (10). With this choice of test functions the integral over
each element Ω
e
vanishes and one is left with integrals on the interfaces between elements and on
the boundary of the domain.
Secondly, we follow the usual idea from finite volume and discontinuous Galerkin methods
of introducing a numerical flux on the interfaces between elements. Consider the interfaces Γ
ee
0
between elements Ω
e
and Ω
e
0
, and on this interface define the unit normal n pointing into Ω
e
0
.
The field variables satisfy the conservation equations (1), and this implies that the flux Fu across
this interface should be continuous. It follows that we can define a numerical flux f
ee
0
such that
f
ee
0
(u
e
, u
e
0
) = F
e
u
e
= F
e
0
u
e
0
. We will discuss the choice of numerical flux in more details in
section 3.3.
Finally we arrive at the following formulation of the wave-based discontinuous Galerkin methods:
X
e,e
0
<e
Z
Γ
ee
0
(v
e
− v
e
0
)
T
f
ee
0
(u
e
, u
e
0
) dΓ +
Z
∂Ω
v
T
Fu dΓ = 0 , ∀v . (12)
The boundary integrals are then modified to implement the different boundary conditions. This
aspect will be discussed in more detail in section 3.4.
3.2. Characteristics
The concept of characteristics plays a central role in the analysis of partial differential equations of
the form (1), see [24], and thus in the construction of numerical fluxes [25]. The basic definitions
and notations are defined in this section to support the discussion of the numerical flux, boundary
conditions and interface conditions in the following sections.
Consider the boundary ∂Ω with unit normal vector n and tangential vector τ . Through a simple
change of variables we can write the governing equations (1) as
iωu + F
∂u
∂n
+ T
∂u
∂τ
= 0 , (13)
where F is the flux matrix defined above and the matrix T = −A
x
n
y
+ A
y
n
x
corresponds to the
flux tangential to boundary.
Copyright
c
2010 John Wiley & Sons, Ltd. Int. J. Numer. Meth. Engng (2010)
Prepared using nmeauth.cls DOI: 10.1002/nme
