Aeroacoustics is a domain where numerical simulation meets great expansion. The min- imization of acoustic pollutions by aircrafts at landing and take o, or more generally by aerospace and terrestrial vehicles, is now an industrial concern, related to more and more severe norms. Dierent approaches coexist under the Computational Aeroacoustics activ- ity. The most widely used methods belong to classical Computational Fluid Dynamics and consist in solving partial dierential equations for the uid, without distinction between the supporting (possibly steady-state) ow and acoustic perturbations. The equations model- ing the uid can be Euler or Navier-Stokes equations, possibly including extended models

clankflaxΜηχανική

22 Φεβ 2014 (πριν από 3 χρόνια και 3 μήνες)

129 εμφανίσεις

Journal of Computational Acoustics,
f
c IMACS
A DISSIPATION-FREE TIME-DOMAIN DISCONTINUOUS GALERKIN METHOD
APPLIED TO THREE-DIMENSIONAL LINEARIZED EULER EQUATIONS
AROUND A STEADY-STATE NON-UNIFORM INVISCID FLOW
MARC BERNACKI
Project-team Caiman,CERMICS and INRIA,
BP93,F06902 Sophia Antipolis Cedex,Marc.Bernacki@sophia.inria.fr
SERGE PIPERNO
Project-team Caiman,CERMICS and INRIA,
BP93,F06902 Sophia Antipolis Cedex,Serge.Piperno@sophia.inria.fr
Received (to be inserted
Revised by Publisher)
We present in this paper a time-domain Discontinuous Galerkin dissipation-free method for the
transient solution of the three-dimensional linearized Euler equations around a steady-state so-
lution.In the general context of a non-uniform supporting ow,we prove,using the well-known
symmetrization of Euler equations,that some aeroacoustic energy satises a balance equation with
source term at the continuous level,and that our numerical framework satises an equivalent bal-
ance equation at the discrete level and is genuinely dissipation-free.In the case of P
1
Lagrange basis
functions and tetrahedral unstructured meshes,a parallel implementation of the method has been
developed,based on message passing and mesh partitioning.Three-dimensional numerical results
conrm the theoretical properties of the method.They include test-cases where Kelvin-Helmholtz
instabilities appear.
Keywords:aeroacoustics;acoustic energy;linearized Euler equations;non-uniformsteady-state ow;
Discontinuous Galerkin method;time domain;energy-conservation.
1.Introduction
Aeroacoustics is a domain where numerical simulation meets great expansion.The min-
imization of acoustic pollutions by aircrafts at landing and take o,or more generally by
aerospace and terrestrial vehicles,is now an industrial concern,related to more and more
severe norms.Dierent approaches coexist under the Computational Aeroacoustics activ-
ity.The most widely used methods belong to classical Computational Fluid Dynamics and
consist in solving partial dierential equations for the uid,without distinction between the
supporting (possibly steady-state) ow and acoustic perturbations
1
.The equations model-
ing the uid can be Euler or Navier-Stokes equations,possibly including extended models
like turbulence,LES techniques,etc
2
.One particular diculty of these approaches is the
dierence in magnitude between the ow and acoustic perturbations,then requiring very
Bernacki,Piperno
accurate { and CPU-consuming { numerical methods.
An alternative has developed recently with approaches consisting in separating the de-
termination of the supporting steady-state ow and in modeling the generation of noise
(for example by providing equivalent acoustic sources),from the propagation of acoustic
perturbations
3;4;5
.For this problem,linearized Euler equations around the supporting ow
are to be solved and provide a good description of the propagation of aeroacoustic pertur-
bations in a smoothly varying heterogeneous and anisotropic medium.This is not exactly
the case of more simple models based on Lighthill analogy
6
or of the third-order equation
of Lilley
7
(a clear description can be found in a more recent reference
8
).The noise source
modeled or derived from the steady-state ow are then dealt with as acoustic source terms
in the linearized Euler equations.
For direct approaches as well as for wave propagation approaches,the construction of
accurate absorbing boundary conditions required for bounding the computational domain
remains a concern.Many solutions have been proposed
10;11;13;14
but the construction of
an general,ecient,parameter-free,easy-to-implement absorbing boundary condition in
time-domain aeroacoustics remains an active research domain
15
.
The work presented here is devoted to the numerical solution of linearized Euler equations
around steady-state discretized ows,obtained using a given Euler solver.The supporting
ow considered is always smooth and subsonic,it can be uniform or fully non-uniform.
Since we intend to consider complex geometries in three space dimensions,we consider
unstructured tetrahedral space discretizations.In this context,we propose a time do-
main Discontinuous Galerkin dissipation-free method based on P
1
Lagrange elements on
tetrahedra.The method is derived from similar methods developed for three-dimensional
time-domain Maxwell equations
16
.We use an element-centered formulation with centered
numerical uxes and an explicit leap-frog time scheme.This kind of method provides a
dissipation-free approximation of propagation equations and allows for the accurate estima-
tion of aeroacoustic energy variation,which is not possible with numerical methods (nite
volumes,discontinuous Galerkin,spectral elements) based on upwind numerical uxes.The
fact that centered numerical uxes in discontinuous Galerkin time-domain methods can lead
to discretization methods inducing no numerical dissipation is quite well known.This was
for example rst numerically shown on cartesian grids
12
,then later for DG methods on
arbitrary unstructured meshes
9;18
.
More precisely,the main results of this paper concern both the linearized Euler equations
at the continuous level,and the numerical DG method we propose.They can be summed
up the following way:
1.for a uniformsupporting ow,at the continuous level (i.e.before space discretization),
some quadratic energy veries a balance equation without source term.This means
energy is conserved (up to boundaries);
2.in this\uniform supporting ow"case,we are able to prove that our DG method
(with leap-frog time-scheme and centered uxes) introduces no dissipation even on
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
unstructured simplicial meshes (some discrete energy is exacly conserved,or simply
non-increasing when absorbing boundary conditions are used);therefore we claimthat
we\control energy variations in the uniform case";
3.accordingly,for a non-uniform supporting ow,at the continuous level (i.e.before
space discretization),we use the well-known symmetrization of nonlinear Euler equa-
tions
17
to derive an aeroacoustic energy which veries some balance equation with
source term.Because of this unsigned source term,aeroacoustic waves can be damped
or excited by the supporting ow.It is responsible for example for Kelvin-Helmholtz
instabilities.These instabilities are due to the model (linearized Euler equations),not
to the numerical method;
4.in the\non uniformsupporting ow"case,we are able to prove that,using an adapted
version of the same DG method on unstructured simplicial meshes,some\discrete"
energy balance equation with source termis also veried.We claim our method is still
non-dissipative.The good point is that we are able to reproduce these instabilities.
The bad point is that we cannot damp themarticially (like methods based on upwind
uxes,which damp instabilities,in an uncontrolled way though);
5.we show nally that,by introducing an articial source term,we are able to\control
these instabilities".
The paper is organized as follows.In Section 2,we present linearized Euler equations
around a steady-state uniformor non-uniformsupporting ow,for which the symmetrization
of Euler equations is introduced to derive some aeroacoustic balance equation.In Section 3,
we present the main elements of the Discontinuous Galerkin time-domain method used,
with the main theoretical results.We give in Section 4 numerical results in two and three
space dimensions in dierent congurations,obtained with MPI-based Fortran parallel im-
plementations of the method.We conclude in Section 5 with possible extensions of this
work.
2.Linearization of Euler equations
We consider in this paper equations for the propagation of acoustic waves through a
steady smooth inviscid ow.Therefore,we linearize the three-dimensional Euler equations
around a given steady ow and only take into account rst-order perturbation terms.For
a perfect inviscid gas,Euler equations read:
@
t
0
B
B
B
B
@

u
v
w
e
1
C
C
C
C
A
+@
x
0
B
B
B
B
@
u
u
2
+p
uv
uw
(e +p)u
1
C
C
C
C
A
+@
y
0
B
B
B
B
@
v
uv
v
2
+p
vw
(e +p)v
1
C
C
C
C
A
+@
z
0
B
B
B
B
@
w
uw
vw
w
2
+p
(e +p)w
1
C
C
C
C
A
= 0;(2.1)
Bernacki,Piperno
where ,~v =
t
(u;v;w),e and p denote respectively the density,the velocity,the volumic
total energy and the pressure,given by the perfect gas law p = ( 1)(e 
1
2
k~vk
2
),where
is a xed constant ( > 1).
2.1.Linearization around a uniform ow
Arst step towards aeroacoustic wave propagation consists in linearizing Euler equations
around any uniform ow dened by some constant physical elds (
0
;~v
0
;p
0
),therefore being
a steady-state solution of Euler equations.Innitely small perturbations (;~v;p) of this
ow verify the following linear hyperbolic system of conservation laws:
@
t
~
W+A
x
@
x
~
W+A
y
@
y
~
W+A
z
@
z
~
W=
~
0;(2.2)
where
~
W=
t
(p c
2
0
;
t

0
c
0
~v;p),c
0
corresponds to the uniform sound speed given by
c
2
0
= p
0
=
0
.The variable
~
Whas been chosen such that the constant matrices A
x
,A
y
,and
A
z
are symmetric.They are given by:
A
x
=
0
B
B
B
B
@
u
0
0 0 0 0
0 u
0
0 0 c
0
0 0 u
0
0 0
0 0 0 u
0
0
0 c
0
0 0 u
0
1
C
C
C
C
A
;A
y
=
0
B
B
B
B
@
v
0
0 0 0 0
0 v
0
0 0 0
0 0 v
0
0 c
0
0 0 0 v
0
0
0 0 c
0
0 v
0
1
C
C
C
C
A
;A
z
=
0
B
B
B
B
@
w
0
0 0 0 0
0 w
0
0 0 0
0 0 w
0
0 0
0 0 0 w
0
c
0
0 0 0 c
0
w
0
1
C
C
C
C
A
In the particular case of the linearization around a uniform ow,one can dene a volumic
aeroacoustic energy E =
1
2
k
~
Wk
2
(a very simple expression in function of
~
W),which veries
the balance equation @
t
E +div
~
F = 0,where the energy ux
~
F is given by
F
s
=
1
2
t
~
WA
s
~
W;s 2 fx;y;zg:
Thus,away from boundaries,the aeroacoustic energy E is exactly conserved.
2.2.Linearization around a non-uniform ow
A more interesting framework consists in linearizing Euler equations around a non-
uniform steady-state solution.In that case,the steady ow is dened by smoothly varying
physical quantities (
0
;~v
0
;p
0
).Linearizing straightforwardly Euler equations (2.1) yields:
@
t
~
W+@
x

A
0
x
~
W

+@
y

A
0
y
~
W

+@
z

A
0
z
~
W

= 0;(2.3)
where
~
Wnow denotes the perturbations of conservative variables (i.e.
~
W
T
= (;
0
~v +
~v
0
;p=( 1) +
0
~v
0
:~v +k~v
0
k
2
=2)) and the space-varying matrices A
0
x
,A
0
y
,and A
0
z
are given in function of ~ = 1,
0
= c
2
0
=~ +k~v
0
k
2
=2,
0
= ( 2)kV
0
k
2
=2 c
2
0
=~ ,and
the canonical basis (~e
x
;~e
y
;~e
z
) of R
3
by
A
0
s
=
0
@
0
t
~e
s
0
~
2
k~v
0
k
2
~e
s
(~v
0
:~e
s
)~v
0
(~v
0
:~e
s
)I
3
~ ~e
s
t
~v
0
+~v
0
t
~e
s
~ ~e
s

0
(~v
0
:~e
s
) 
0
t
~e
s
 ~ (~v
0
:~e
s
)
t
~v
0
(~v
0
:~e
s
)
1
A
;s 2 fx;y;zg:(2.4)
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
In this equation,the matrices A
0
x
,A
0
y
,and A
0
z
are not symmetric anymore,and it is
very dicult to deduce any aeroacoustic energy balance equation.Therefore,we con-
sider other acoustic variables,derived from the quite classical symmetrization of Euler
equations.Assuming the ow is smooth enough,the change of variables (;~v;e)!
(
e~
p
+ +1 ln

p



;
~
p
~v;
~
p
) transforms Euler equations (2.1) into a symmetric sys-
tem of conservation laws (i.e.Jacobians of uxes are symmetric matrices).Accordingly,
the linearization of these symmetrized Euler equations leads to more complex aeroacoustic
equations for perturbations of the new variables,which can be written as
A
0
0
@
t
~
V+@
x

~
A
0
x
~
V

+@
y

~
A
0
y
~
V

+@
z

~
A
0
z
~
V

= 0;(2.5)
where
~
V is given in function of variables
~
Wby
~
V = A
0
0
1
~
Wand
A
0
0
=

0
~
0
B
@
1
t
~v
0

0
c
2
0
=
~v
0
c
2
0

I
3
+~v
0
t
~v
0

0
~v
0

0
c
2
0
= 
0
t
~v
0

2
0
c
4
0
= ( ~ )
1
C
A
;
~
A
0
s
= A
0
s
A
0
0
;s 2 fx;y;zg:(2.6)
A
0
0
clearly is symmetric and it can be proved that it is denite positive (and then not
singular).Eq.2.5 can also be obtained simply by replacing
~
Wby A
0
0
~
V in Eq.2.3 (and by
noting that @
t
A
0
0
= 0).Finally,the reader can also check that the symmetric matrices
~
A
s
(s 2 fx;y;zg) are given by
~
A
0
s
= (~v
0
:~e
s
) A
0
0
+
p
0
~
0
@
0
t
~e
s
(~v
0
:~e
s
)
~e
s
~e
s
t
~v
0
+~v
0
t
~e
s
(~v
0
:~e
s
)~v
0
+
0
~e
s
(~v
0
:~e
s
) (~v
0
:~e
s
)
t
~v
0
+
0
t
~e
s
2
0
(~v
0
:~e
s
)
1
A
:
Then,the volumic aeroacoustic energy E dened by E =
1
2
t
~
WA
0
0
1
~
W
1
2
t
~
VA
0
0
~
V veries
the following balance equation with source term:
@
t
E +div
~
F = S;with
(
F
s
=
t
~
V
~
A
0
s
~
V;s 2 fx;y;zg:
S = 
1
2
t
~
V
h
@
x
(
~
A
0
x
) +@
y
(
~
A
0
y
) +@
z
(
~
A
0
z
)
i
~
V:
(2.7)
Thus the aeroacoustic energy is not conserved and the variations in the steady ow con-
sidered can damp or amplify aeroacoustic waves,unless the source term vanishes (which is
the case for a uniform ow for example).In the sequel,we shall mainly discretize the con-
servative form (2.3),but we shall need the equivalent symmetric form (2.5) for discussions
concerning energy conservation and stability.
3.A discontinuous Galerkin time-domain (DGTD) method
Discontinuous Galerkin methods have been widely used with success for the numerical
simulation of acoustic or electromagnetic wave propagation in the time domain
16;18
.The
very same type of methods can be used for the problems considered here,i.e.the propagation
Bernacki,Piperno
of aeroacoustic waves through a non-uniform ow
19;20;21;22
.In this section,we present the
DGTD method we use for the model equations (2.3).We recall the numerical properties
of the space discretization.Then we introduce the leap-frog time scheme and give some
details on properties related to energy conservation and stability.
In the whole paper,we assume we dispose of a partition of a polyhedral domain

(whose boundary @
is the union of physical boundaries of objects @

phys
and of the far
eld articial boundary @

1
).
is partitioned into a nite number of polyhedra (each
one having a nite number of faces).For each polyhedron T
i
,called"control volume"or
"cell",V
i
denotes its volume.We call face between two control volumes their intersection,
whenever it is a polyhedral surface.The union of all faces F is partitioned into internal faces
F
int
= F=@
,physical faces F
phys
= F
T
@

phys
and absorbing faces F
abs
= F
T
@

1
.For
each internal face a
ik
= T
i
T
T
k
,we denote by S
ik
the measure of a
ik
and by ~n
ik
the unitary
normal,oriented from T
i
towards T
k
.The same denitions are extended to boundary faces,
the index k corresponding to a ctitious cell outside the domain.Finally,we denote by V
i
the set of indices of the control volumes neighboring a given control volume T
i
(having a
face in common).We also dene the perimeter P
i
of T
i
by P
i
=
P
k2V
i
S
ik
.We recall the
following geometrical property for all control volumes:
P
k2V
i
S
ik
~n
ik
= 0.
Following the general principle of discontinuous Galerkin nite element methods,the
unknown eld inside each control volume is seeked for as a linear combination of local
basis vector elds ~'
ij
;1  j  d
i
(generating the local space P
i
) and the approximate
eld is allowed to be fully discontinuous across element boundaries.Thus,a numerical ux
function has to be dened to approximate uxes at control volumes interfaces,where the
approximate solution is discontinuous.
This context is quite general.Actual implementations of the method have been consid-
ered only on tetrahedral meshes,where control volumes are the tetrahedra themselves.We
shall only consider constant (P
0
) or linear (P
1
) approximations inside tetrahedra.
3.1.Time and space discretizations
We only consider here the most general case of aeroacoustic wave propagation in a non-
uniform steady ow.Also,in order to limit the amount of computations,we restrict our
study to piecewise constant matrices A
0
s
(s 2 fx;y;zg) given in Eq.(2.4).For each control
volume T
i
,for s 2 fx;y;zg,we denote by A
i
s
an approximate for the average value of A
0
s
over T
i
.Dot-multiplying Eq.(2.3) by any given vector eld ~',integrating over T
i
and
integrating by parts yields
Z
T
i
~'
@
~
W
@t
=
Z
T
i
0
@
X
s2fx;y;zg
t
@
s
~'A
0
s
1
A ~
W
Z
@T
i
~'
0
@
X
s2fx;y;zg
n
s
A
0
s
~
W
1
A
:(3.8)
Inside volume integrals over T
i
,we replace the eld
~
W by the approximate eld
~
W
i
and
the matrices A
0
s
by their respective average values A
i
s
.For boundary integrals over @T
i
,
~
W
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
is discontinuous,and we dene totally centered numerical uxes,i.e.:
(
8i;8k 2 V
i
;
h
(n
ik
x
A
0
x
+n
ik
y
A
0
y
+n
ik
z
A
0
z
)
~
W
i
ja
ik
'
1
2

P
i
ik
~
W
i
+P
k
ik
~
W
k

;
with P
i
ik
= n
ik
x
A
i
x
+n
ik
y
A
i
y
+n
ik
z
A
i
z
;P
k
ik
= n
ik
x
A
k
x
+n
ik
y
A
k
y
+n
ik
z
A
k
z
:
(3.9)
Concerning the time discretization,we use a three-level leap-frog scheme.The unknowns
~
W
i
are approximated at integer time-stations t
n
= nt.Assuming we dispose of
~
W
n1
i
and
~
W
n
i
,the unknowns
~
W
n+1
i
are seeked for in P
i
such that,8~'2 P
i
,
Z
T
i
~'
~
W
n+1
i

~
W
n1
i
2t
=
Z
T
i
X
s2fx;y;zg
t
@
s
~'A
i
s
~
W
n
i

X
k2V
i
Z
a
ik
~'
P
i
ik
~
W
n
i
+P
k
ik
~
W
n
k
2
:(3.10)
Again,the time scheme above is almost explicit.Each time step only requires the inversion
of local symmetric positive denite matrices of size (d
i
d
i
).In the particular case where
P
i
is a complete linear (P
1
) representation,these 20  20 matrices are indeed made of 5
4 4 diagonal blocks (which are equal to the classical P
1
mass matrix).
3.2.Boundary conditions
Boundary conditions are dealt with in a weak sense.For the physical boundary,we
consider only a slip condition,which is set on both the steady ow and the acoustic per-
turbations.This means that we assume that for any slip boundary face a
ik
belonging to
the control volume T
i
,the steady solution of Euler equations veries a slip condition at
the discrete level,i.e.~n
ik
 ~v
0
i
= 0.For the acoustic perturbations,we use a mirror c-
titious state
~
W
k
in the computation of the boundary ux given in Eq.(3.10).We take

k
= 
i
,p
k
= p
i
,and ~v
k
= ~v
i
2(~n
ik
 ~v
i
)~n
ik
(which implies (~v
k
~v
i
) ~n
ik
= 0 and
~v
k
:~n
ik
= ~v
i
:~n
ik
).
For an absorbing boundary face a
ik
,upwinding is used to select outgoing waves only.
Before discretization in time,classical upwinding leads to a boundary ux F
ik
given by
F
ik
=

P
i
ik

+
~
W
i
,where for any diagonalizable matrix Q = S
1
DS with D diagonal,Q
+
=
(Q+jQj)=2 and terms of the diagonal matrix jDj are the moduli of the eigenvalues.This
general idea leads to P
k
ik
~
W
k
= jP
i
ik
j
~
W
i
.However,for this intuitive numerical ux,it is
very dicult to prove that the resulting time-scheme is stable and that energy is actually
sent in the exterior domain.We then consider the numerical ux based on the following
ctitious state:P
k
ik
~
W
n
k
=
p
A
i
0



p
A
i
0
1
P
i
ik
p
A
i
0



p
A
i
0
1
~
W
n1
i
+
~
W
n+1
i
2
,where
p
A
i
0
is the positive
square root of the symmetric denite positive matrix A
i
0
.Indeed,this expression derives
from the intuitive upwind ux for the symmetrized equations (2.5).It leads to time-scheme
which is locally implicit near absorbing boundaries (i.e.independent linear systems are to
be solved inside elements having at least one absorbing face,at each time step).It leads
to a globally second-order time-accurate scheme.A less accurate explicit version is also
available
22
.It takes the form P
k
ik
~
W
n
k
=
p
A
i
0



p
A
i
0
1
P
i
ik
p
A
i
0



p
A
i
0
1
~
W
n1
i
.
Bernacki,Piperno
3.3.Energy balance and stability
In order to investigate stability,we dene a discrete aeroacoustic energy F
n
by:
F
n
=
1
2
X
i
Z
T
i
t
~
W
n
i
A
i
0
1
~
W
n
i
+
t
~
W
n+1
i
A
i
0
1
~
W
n1
i
+
t
4
X
a
ik
2F
abs
Z
a
ik
t

A
i
0
1
~
W
n1
i

M
ik

A
i
0
1

~
W
n1
i
+
~
W
n+1
i

;
with M
ik
=
p
A
i
0



p
A
i
0
1
P
i
ik
p
A
i
0



p
A
i
0

p
A
i
0



p
A
i
0
1
~
P
i
ik
p
A
i
0
1



p
A
i
0
.One can show that the
matrices M
ik
are symmetric and positive.We give in Annex 1 the expression of the variation
through one time step of the aeroacoustic energy (i.e.F
n+1
F
n
):
F
n+1
F
n
= 
t
2
X
a
ik
2F
int
Z
a
ik
t
~
V
n
i
(
~
P
k
ik

~
P
i
ik
)
~
V
n+1
k
+
t
~
V
n+1
i
(
~
P
k
ik

~
P
i
ik
)
~
V
n
k

t
4
X
a
ik
2F
abs
Z
a
ik
t

~
V
n1
i
+
~
V
n+1
i

M
ik

~
V
n1
i
+
~
V
n+1
i

:(3.11)
where we have used the auxiliary variables
~
V
n
i
 A
i
0
1
~
W
n
i
;8i;8n.The rst term is a
discrete version of the source term appearing in Eq.(2.7).The second term is negative and
shows that our absorbing boundary conditions actually absorbs energy.This results also
shows that the slip boundary condition has no in uence on the global energy balance.
In order to prove stability,we show in Annex 2 that F
n
is a quadratic positive denite
formof numerical unknowns (
~
W
n1
i
;
~
W
n
i
) under some CFL-like sucient stability condition
on the time-step t:
8i;8k 2 V
i
;t (2
i

i
+
ik

ik
) <
2V
i
P
i
;(3.12)
where 
i
and 
ik
are dimensionless regularity coecients depending of basis functions and
element aspect ratio,
i
= ju
i
0
j +jv
i
0
j +jw
i
0
j +3c
i
0
,and 
ik
= j~v
i
0
~n
ik
j +c
i
0
for a boundary face
and 
2
ik
= sup


j~v
i
0
 ~n
ik
j +c
i
0
)

2


A
k
0
A
i
0
1

;

j~v
k
0
 ~n
ik
j +c
k
0

2


A
i
0
A
k
0
1

for an internal
face ( here denotes the spectral radius of a matrix).
In the case of a uniform ow,we have P
i
ik
= P
k
ik
.Thus the aeroacoustic energy is
non-increasing (and exactly conserved if no absorbing boundary is present,which shows
the scheme is genuinely non-diusive) and the scheme is stable under a CFL-type stability
condition depending on the size of elements and sup
i
(k~v
i
k +c
0
i
).
4.Numerical results
We dispose of a three-dimensional parallel implementation of the DGTD method pre-
sented in the previous section.Any subsonic steady ow can be considered,even with strong
ow gradients.However,the ow,given as the output of a non-linear Euler equations solver,
has to be post-processed:average of the ow over tetrahedra must be computed and the
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
non-slip condition must be enforced on physical boundaries.We present in this section test-
cases in two and three space dimensions,in order to validate the method on benchmark
problems,test the method on complex ows and congurations,and nally evaluate the
performance of the parallel Fortran 77 implementation,based on the MPICH implementa-
tion of MPI.Parallel computations were performed on a 16 node cluster (2GHz-Pentium4
1Gb-RDRAM memory biprocessor each).In this section,tables give performance results
for 64 bit arithmetic computations:N
p
is the number of processes for the parallel execution,
REAL denotes the total (wall clock) simulation time and CPU denotes the corresponding
total CPU time taken as the maximum of the per process values.Finally,% CPU denotes
the ratio of the total CPU time to the total wall clock time.This ratio clearly allows an
evaluation of the CPU utilization and yields a metric for parallel eciency.
4.1.Acoustic perturbations in a three-dimensional uniform ow
This quite classical test-case
23;24
is aimed at validating the DGTD method proposed by
mixing acoustic,entropic and vorticity perturbations,in a uniform horizontal ow (with
Mach number 0.5).In that case,the exact solution is well known.The geometry is a unit
cube with absorbing boundaries.Two unstructured tetrahedral meshes have been used,
whose characteristics are given in Table 1 below.Plane cuts in the numerical solution
Table 1.Tam's test-case.Characteristics of meshes
Mesh#vertices#tetrahedra#absorbing faces
M1 68,921 384,000 19,200
M2 531,441 3,072,000 76,800
are shown on Fig.1.The result agree with the exact solution.The absorbing boundary
condition is not perfect,but the L2-norm of the error at time t = 50 is limited to 2% for
the coarse mesh M1 and around 1% for the ner mesh M2.Articial re ections at the
absorbing boundaries are limited to 1% (of the amplitude at the corresponding time) for
mesh M1 and around 0;5% for M2.The time evolution of the energy inside the domain is
given on Fig.2.As expected,since the steady ow is uniform,the total acoustic energy is
rst exactly conserved,and then decreasing,once perturbations have reached the absorbing
boundaries (at t = 15).Computation times are given in Table 2.Parallel eciency reaches
a satisfying level and total CPU time are really acceptable.
4.2.Acoustic perturbations in a three-dimensional subsonic Couette ow
This second test-case is more complex and is classically used to evaluate the merits of a
numerical method for non-uniform shear ows
25
.The steady ow considered is horizontal,
with ~v
0
=
t
(0:9z;0;0).At the beginning of the simulation,the same (acoustic,entropic,
and vorticity) perturbations are mixed.The domain considered is parallelepipedic with a
larger dimension in the direction of the supporting ow.Although the geometry recalls a
Bernacki,Piperno
Fig.1.Tam's test-case. contour lines at times t = 10:7,25,39:3,and 50.
0.5
1
1.5
2
2.5
3
3.5
0
5
10
15
20
25
30
35
40
45
50
55
F
time
Fig.2.Tam's test-case.Time evolution of acoustic energy.
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
Table 2.Tam's test-case.Parallel CPU times and eciency.
Mesh N
p
CPU REAL % CPU S(N
p
)
M1 1 33381 s 33447 s 100% 1
- 4 8167 s 8412 s 97% 4
- 8 4286 s 4465 s 96% 7.5
- 16 2195 s 2381 s 92% 15.2
M2 8 34012 s 35651 s 95% 1
- 16 17203 s 18760 s 92% 1.98
waveguide,the aim of this test-case is to validate the accurate propagation and convection
of the waves by the shear ow,and the behavior of the absorbing boundary condition,which
is set on all boundaries of the domain.Two unstructured tetrahedral meshes have been
generated.Their characteristics are given in Table 3.Plane cuts and volumic contours for
Table 3.Couette ow test-case.Characteristics of meshes.
Mesh#vertices#tetrahedra#absorbing faces
M1 73,629 405,600 23,504
M2 522,801 3,000,000 90,000
 in the numerical solution are shown on Fig.3.The accuracy of the numerical results
is very similar to the one observed for a uniform ow:in the numerical results,articial
re ections fromabsorbing boundaries have a relative amplitude of 1:5% for the coarse mesh
M1 (respectively 1% for the ner mesh M2) compared to maximal amplitude at the same
time in the whole domain.In this test-case,the initial acoustic perturbation propagates
and gets quickly out of the domain,while the entropic perturbation is only convected by
the supporting ow (thus the bottom part of this perturbation travels slowly with the
supporting ow,whose velocity vanishes in the plane z = 0).
The time evolution of the global acoustic energy is given on Fig.4.As expected,the en-
ergy is not conserved anymore because the non-uniform ow provides energy to the acoustic
perturbations.The energy decreases when perturbations reach the absorbing boundaries.
A part of the energy corresponding to the bottompart of the entropic perturbation remains
in the domain.Computation times are given in Table 4 and reveal similar parallel eciency
and execution times.
4.3.Two-dimensional Kelvin-Helmholtz instabilities
It is well-known that particular steady-state solutions of Euler equations can lead to
unstable solutions for the corresponding linearized Euler equations.These instabilities are
known as Kelvin-Helmholtz instabilities and are proved to appear for example for shear
ows.These instabilities are present in the linearized equations.In the original Euler
equations,some limiting nonlinear term must be playing the role of a limiter for the overall
Bernacki,Piperno
Fig.3.Couette ow test-case. contour lines (left:contours in plane cuts;right:volumic contours) at
times t = 0,0:5,and 5:75.
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
0.08
0
1
2
3
4
5
Energy
Time (s)
Fig.4.Couette ow test-case.Time evolution of the acoustic energy.
Table 4.Couette ow test-case.Parallel CPU times and eciency.
Mesh N
p
CPU REAL % CPU S(N
p
)
M1 4 12830 s 13008 s 99% 1
- 8 6428 s 6590 s 98% 2
- 16 3230 s 3373 s 96% 3.9
M2 8 46359 s 47743 s 97% 1
- 16 23412 24959 s 94% 1.98
Bernacki,Piperno
acoustic energy in the ow.When linearized Euler equations are solved using numerical
methods with numerical dissipation,these instabilities might be eliminated.It might be an
positive feature in some cases.However,there is no control on the numerical dissipation and
the accuracy of the solution obtained.In our case,the numerical method we have proposed
is genuinely non-dissipative.Therefore,Kelvin-Helmholtz instabilities should appear easily,
and more easily as transverse gradients in the steady ow increase.This is the case for
example in two space dimensions with ~v
0
= 0:25 tanh [151:51(y 0:25) +0:5] ~e
x
and 
0
= 1,
p
0
= 1= .In this ow,u
0
goes through an in ection point at y = 0:25 (and k~v
0
k = 0:5
on the in ection point).In this particular ow,the acoustic energy is amplied along the
in ection line and at the same time,perturbations are convected along this line,which
leads to an instability.The behavior is completely dierent from the one obtained with a
shear ow of the form ~v
0
= (y +0:25)~e
x
,where no instability occurs.We have introduced
a periodic in time,Gaussian in space source term to excite possible unstable modes and
numerical results obtained are compared on Fig.5.The instability appears quite quickly
for the shear ow with an in ection point.Dierent means are under current investigation
to control this kind of instability.
Fig.5.Contours of k~vk at t=0.257 for the linear (left) and in ection (right) shear ows with source term
4.4.Two-dimensional instabilities past a NACA airfoil prole
This test-case is aimed at evaluating the occurrence of Kevin-Helmholtz-type instabilities
for steady-state ows with sharp gradients.We have considered a test-case past a NACA-
type airfoil prole (courtesy of ONERA
19
).The steady-state solution at M
1
= 0:5 obtained
by our nite-volume non-linear Euler solver on the 33046-vertex 65580-triangle unstructured
mesh is shown on Fig.6.A Gaussian in space and periodic in time (T = 1ms) source term
has been set up for the aeroacoustic problem,where a slip boundary condition is used on
the prole,and our absorbing boundary conditions in the far eld boundary.For this test-
case,we have observed an instability,i.e.the energy source term has a positive sign and
some aeroacoustic blow up occurs.A detail on contours of k~vk are shown on Fig.7.Some
numerical treatments are possible to overcome these diculties and get rid of instabilities.
One of them
4
consists in adding a source term which leads to energy dissipation.As we
are not specialists of physics,we do not discuss here the validity of such a model,but we
can observe that the instability disappears after treatment.Contours of k~vk are shown on
the right part of Fig.8 and can be compared to classical acoustic scattering (uniform still
ow).The comparison of contours near the trailing edge for acoustics and aeroacoustics
is shown on Fig.9,where the shear ow at the trailing edge has a visible in uence on the
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
Mach, min = 0.0231943, max = 0.889346
Fig.6.Mach number contours for the steady-state ow past the airfoil prole at M
1
= 0:5
Fig.7.Contours of k~vk at t = 1:57ms for the original aeroacoustic model.
Bernacki,Piperno
propagation of acoustics (and generation of vortices).
Fig.8.Contours of k~vk for the stabilized aeroacoustic model (left) and classical acoustics (right).
Fig.9.Tailing edge zoom { k~vk contours for stabilized aeroacoustics (left) and acoustics (right).
4.5.Three-dimensional Aeroacoustic wave propagation past a sphere
We consider the case of a steady three-dimensional ow past a sphere and the propa-
gation of waves emitted by an acoustic source (Gaussian prole,periodic source term with
period T = 0:2s.The subsonic ow with M
1
= 0:5 is obtained by a three-dimensional
MUSCL-extended nite-volume solver of three-dimensional Euler equations on tetrahedral
grids.The steady-state solution presents Mach numbers between 0:002 and 0:8 and con-
tours for 
0
are presented on Fig.10.Details on the computational tetrahedral mesh are
given on Table 5.The unit sphere is centered inside the 10m-side cubic computational
domain.A slip boundary condition is set on the sphere,while the absorbing boundary
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
Figure 10:Aeroacoustic propagation past a sphere { steady-state density 
0
condition is set on far eld boundaries.The aeroacoustic simulation was performed on 16
processors,with computational times reported on Table 5.We have shown on Fig.11 the
Table 5:Aeroacoustic propagation past a sphere { computational mesh and times
#vertices#elements#abs.faces#slip faces N
p
CPU REAL % CPU
324,471 1,870,288 10,092 46,080 16 50h4mn 51h18mn 97,6%
resulting propagated waves for this aeroacoustic case in a steady ow past a sphere (after
10 and 25 periods).The numerical results (scattering and re ection patterns) as well as
the computational times are encouraging.One can notice that the periodic regime is not
reached after 10 periods.These results lead us to consider in the near future more complex
test-cases and more realistic congurations.
5.Conclusion
In this paper,we have presented a time domain Discontinuous Galerkin dissipation-
free method for the transient solution of the three-dimensional linearized Euler equations
around a steady-state solution,based on P
1
Lagrange elements on tetrahedra.In the more
general context of a non-uniform supporting ow,we have proved,using the well-known
symmetrization of Euler equations,that the aeroacoustic energy satises a balance equa-
tion with source term at the continuous level,and that our numerical method satises
an equivalent balance equation at the discrete level.This shows that the method is gen-
uinely dissipation-free,which is conrmed by numerical results,including test-cases where
Kelvin-Helmholtz instabilities appear.
Further works can concern many dierent aspects.On the modeling side,it is possible to
Bernacki,Piperno
Figure 11:Aeroacoustic propagation past a sphere {  after 10T (left) and 25T (right)
design models for dealing with natural Kevin-Helmholtz instability,for example by adding
some source terms.This has to be done in cooperation with physicists.Anyway,the
numerical framework proposed here provides a valuable tool for investigating this kind of
instability in complex ows and geometries.On the numerical side,the overall accuracy
could be enhanced either by considering more-than-linear basis functions (P
k
Lagrange
elements with k > 1) or by dealing with a more accurate description of the supporting
ow (currently,it is only P
0
).Higher-order accuracy in absorbing boundary conditions and
time-scheme should also be seeked for.
Annex 1
We shall prove the result given in Eq.(3.11).Since A
i
0
is symmetric,we have
F
n+1
F
n
=
1
2
X
i
Z
T
i
t
~
W
n
i
A
i
0
1

~
W
n+2
i

~
W
n
i

+
t
~
W
n+1
i
A
i
0
1

~
W
n+1
i

~
W
n1
i

+
t
4
X
a
ik
2F
abs
Z
a
ik

t

A
i
0
1
~
W
n
i

M
ik
A
i
0
1

~
W
n+2
i
+
~
W
n
i


t

A
i
0
1
~
W
n1
i

M
ik
A
i
0
1

~
W
n+1
i
+
~
W
n1
i


:
The variational formof the scheme (3.10) allows to consider any eld in P
i
,thus for example
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
the elds
~
W
n1
i
or
~
W
n
i
themselves.We obtain that F
n+1
= F
n
+t(T
1

1
2
T
2

1
2
T
3
) with
T
1
=
X
i
Z
T
i
X
s2fx;y;zg
t
@
s

~
V
n
i

~
A
i
s
~
V
n+1
i
+
t
@
s

~
V
n+1
i

~
A
i
s
~
V
n
i
T
2
=
X
a
ik
2F
int
Z
a
ik
t
~
V
n
i

~
P
i
ik
~
V
n+1
i
+
~
P
k
ik
~
V
n+1
k

+
t
~
V
n
k

~
P
k
ki
~
V
n+1
k
+
~
P
i
ki
~
V
n+1
i

+
X
a
ik
2F
int
Z
a
ik
t
~
V
n+1
i

~
P
i
ik
~
V
n
i
+
~
P
k
ik
~
V
n
k

+
t
~
V
n+1
k

~
P
k
ki
~
V
n
k
+
~
P
i
ki
~
V
n
i

T
3
=
X
a
ik
2F
slip
Z
a
ik
t
~
V
n
i

~
P
i
ik
~
V
n+1
i
+
~
P
i
ik
H
ik
~
V
n+1
i

+
t
~
V
n+1
i

~
P
i
ik
~
V
n
i
+
~
P
i
ik
H
ik
~
V
n
i

+
X
a
ik
2F
abs
Z
a
ik
"
t
~
V
n
i

~
P
i
ik
~
V
n+1
i
+M
ik
~
V
n
i
+
~
V
n+2
i
2
!
+
t
~
V
n+1
i

~
P
i
ik
~
V
n
i
+M
ik
~
V
n1
i
+
~
V
n+1
i
2
!#

1
2
X
a
ik
2F
abs
Z
a
ik
t
~
V
n
i
M
ik

~
V
n
i
+
~
V
n+2
i


t
~
V
n1
i
M
ik

~
V
n1
i
+
~
V
n+1
i

:
Since
~
P is symmetric,and since
~
P
i
ik
= 
~
P
i
ki
and
~
P
k
ik
= 
~
P
k
ki
,we obtain that
T
2
= 2
X
a
ik
2F
int
Z
a
ik
t
~
V
n
i
~
P
i
ik
~
V
n+1
i
+
t
~
V
n
k
~
P
k
ki
~
V
n+1
k
+
X
a
ik
2F
int
Z
a
ik
t
~
V
n
i

~
P
k
ik

~
P
i
ik

~
V
n+1
k
+
t
~
V
n
k

~
P
k
ik

~
P
i
ik

~
V
n+1
i
;and
T
2
+T
3
= 2
X
i
X
k2V
i
Z
a
ik
t
~
V
n
i
~
P
i
ik
~
V
n+1
i
+
X
a
ik
2F
slip
Z
a
ik
t
~
V
n
i
~
P
i
ik
H
ik
~
V
n+1
i
+
t
~
V
n+1
i
~
P
i
ik
H
ik
~
V
n
i
+
X
a
ik
2F
int
Z
a
ik
t
~
V
n
i

~
P
k
ik

~
P
i
ik

~
V
n+1
k
+
t
~
V
n
k

~
P
k
ik

~
P
i
ik

~
V
n+1
i
+
1
2
X
a
ik
2F
abs
Z
a
ik
t

~
V
n1
i
+
~
V
n+1
i

M
ik

~
V
n1
i
+
~
V
n+1
i

:
In the expressions above,the matrix H
ik
denotes the matrix used for slip boundary con-
ditions,in the sense that the ctitious state
~
W
n
k
is given by
~
W
n
k
= H
ik
~
W
n
i
(H
ik
does not
change the density and energy perturbations,and operates an orthogonal symmetry on the
velocity perturbation).One can also show that H
2
ik
= I and that H
ik
commutes with A
i
0
.
Bernacki,Piperno
We nally get that F
n+1
= F
n
+t(T
4
+T
5
+T
6
) with
T
4
=
X
i
2
4
Z
T
i
X
s2fx;y;zg

t
@
s

~
V
n
i

~
A
i
s
~
V
n+1
i
+
t
@
s

~
V
n+1
i

~
A
i
s
~
V
n
i


X
k2V
i
Z
a
ik
t
~
V
n
i
~
P
i
ik
~
V
n+1
i
3
5
T
5
= 
1
2
X
a
ik
2F
ref
Z
a
ik
t
~
V
n
i

~
P
i
ik
H
ik
+
t

~
P
i
ik
H
ik


~
V
n+1
i
T
6
= 
1
4
X
a
ik
2F
abs
Z
a
ik
t

~
V
n1
i
+
~
V
n+1
i

M
ik

~
V
n1
i
+
~
V
n+1
i


1
2
X
a
ik
2F
int
Z
a
ik
t
~
V
n
i

~
P
k
ik

~
P
i
ik

~
V
n+1
k
+
t
~
V
n
k

~
P
k
ik

~
P
i
ik

~
V
n+1
i
Since the matrices
~
A
s
are symmetric and constant over control volumes,we get by integrat-
ing by parts that T
4
= 0.Also,since the matrix
~
P
i
ik
H
ik
is skew-symmetric,we get T
5
= 0.
Using the expression of M
ik
and using the
~
Wvariables,we get the result announced.
Annex 2
We shall prove that F
n
is a positive denite quadratic form of unknowns
~
W
n1
i
and
~
W
n
i
under the condition given in Eq.(3.12).Let us introduce the following denitions.
Denition 1 8
~
X 2 P
i
,we denote by k
~
Xk
T
i
= (
R
T
i
k
~
Xk
2
)
1=2
the L
2
-norm of
~
X over T
i
.
For any interface a
ik
2 F,we also dene k
~
Xk
a
ik
= (
R
a
ik
k
~
Xk
2
)
1=2
.Finally,the matrix A
i
0
being positive denite,we dene a second norm by kj
~
Xkj
T
i
= k
p
A
i
0
1
~
Xk
T
i
.
Denition 2 We assume some regularity for the basis functions ~'
ij
;1  j  d
i
.More
precisely,we assume that,for each control volume T
i
,there exist dimensionless positive
constants 
i
and 
ik
(for each k 2 V
i
) such that
8
~
X2 P
i
;
(
8s 2 fx;y;zg;k@
s
~
Xk
T
i


i
P
i
V
i
k
~
Xk
T
i
;
8a
ik
2 F;k
~
Xk
2
a
ik


ik
k~n
ik
k
V
i
k
~
Xk
2
T
i
:
(5.13)
Denition 3 We choose here some notations.We dene
~
Z
n
i
=
p
A
i
0
1
~
W
n
i
.We introduce

i
=
P
s


A
i
s

(and we recall 

A
i
s

= j~e
s
:~v
i
0
j + c
i
0
).For boundary faces a
ik
,we recall
that some neighboring ctitious control volume T
k
is imagined,and we take by convention:
k
~
W
n
k
k
T
k
= k
~
W
n
i
k
T
i
,kj
~
W
n
k
kj
T
k
= kj
~
W
n
i
kj
T
i
,
ki
= 
ik
,V
k
= V
i
.Finally,we dene 
ik
for
any face a
ik
by:
8
>
<
>
:
a
ik
2 F
int
;
2
ik
= sup


j~v
i
0
 ~n
ik
j +c
i
0

2


A
k
0
A
i
0
1

;

j~v
k
0
 ~n
ik
j +c
k
0

2


A
i
0
A
k
0
1

a
ik
2 F
slip
;
ik
= 

P
i
ik

= j~v
i
0
 ~n
ik
j +c
i
0
a
ik
2 F
abs
;
ik
= 0;
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
In order to prove that F
n
is a positive denite quadratic form of unknowns
~
W
n1
i
and
~
W
n
i
under the condition (3.12),we rst show the following result:
F
n
i

1
2

kj
~
W
n
i
kj
2
T
i
+kj
~
W
n1
i
kj
2
T
i



i
P
i

i
t
V
i
kj
~
W
n
i
kj
T
i
kj
~
W
n1
i
kj
T
i

t
4
X
k2V
i


ik

ik
k~n
ik
k
V
i
kj
~
W
n1
i
kj
2
T
i
+

ki

ki
k~n
ki
k
V
k
kj
~
W
n
k
kj
2
T
k

:(5.14)
We can show easily that the above lemma gives the result.Using the above denitions and
assumptions,we have F
n
i

P
k2V
i
k~n
ik
k(X
ik


ki

ki
t
4V
k
kj
~
W
n
k
kj
2
T
k
) with
X
ik
=
1
2P
i
kj
~
W
n
i
kj
2
T
i
+

1
2P
i


ik

ik
t
4V
i

kj
~
W
n1
i
kj
2
T
i


i

i
t
V
i
kj
~
W
n
i
kj
T
i
kj
~
W
n1
i
kj
T
i


1
2P
i


i

i
t
2V
i

kj
~
W
n
i
kj
2
T
i
+

1
2P
i


ik

ik
t
4V
i


i

i
t
2V
i

kj
~
W
n1
i
kj
2
T
i
Finally,F
n

X
a
ik
2F
int
k~n
ik
kY
ik
+
X
a
ik
2F
slip
S
F
abs
k~n
ik
kZ
ik
with
Y
ik
=

1
2P
i


ik

ik
t
4V
i


i

i
t
2V
i


kj
~
W
n
i
kj
2
T
i
+kj
~
W
n1
i
kj
2
T
i

+

1
2P
k


ki

ki
t
4V
k


k

k
t
2V
k


kj
~
W
n
k
kj
2
T
k
+kj
~
W
n1
k
kj
2
T
k

;
Z
ik
=

1
2P
i


ik

ik
t
4V
i


i

i
t
2V
i


kj
~
W
n
i
kj
2
T
i
+kj
~
W
n1
i
kj
2
T
i

:
Therefore,the energy F
n
is a positive denite quadratic form of unknowns under the su-
cient CFL-type stability condition proposed in Eq.(3.12).
End of the proof.We now give the proof for the preliminary result (5.14).We rst have
F
n
i
=
1
2

kj
~
W
n
i
kj
2
T
i
+kj
~
W
n1
i
kj
2
T
i

+
t
2
(X
n
i
+Y
n
i
) with
X
n
i
=
1
2
X
a
ik
2F
abs
\@T
i
Z
a
ik
t

A
i
0
1
~
W
n1
i

M
ik
A
i
0
1

~
W
n1
i
+
~
W
n+1
i

=
1
2
X
a
ik
2F
abs
\@T
i
Z
a
ik
t
~
V
n1
i
M
ik

~
V
n1
i
+
~
V
n+1
i

Y
n
i
=
1
t
Z
T
i
t
~
W
n1
i
A
i
0
1

~
W
n+1
i

~
W
n1
i

= 2
Z
T
i
X
s2fx;y;zg
t
@
s

~
V
n1
i

~
A
i
s
~
V
n
i

X
k2V
i
Z
a
ik
t
~
V
n1
i

~
P
i
ik
~
V
n
i
+
~
P
k
ik
~
V
n
k

=
Z
T
i
X
s2fx;y;zg

t
@
s

~
V
n1
i

~
A
i
s
~
V
n
i

t
~
V
n1
i
~
A
i
s
@
s
~
V
n
i


X
k2V
i
Z
a
ik
t
~
V
n1
i
~
P
k
ik
~
V
n
k
Bernacki,Piperno
Using H
ik
and eliminating terms corresponding to absorbing boundaries leads to
X
n
i
+Y
n
i
=
Z
T
i
X
s2fx;y;zg

t
@
s

~
V
n1
i

~
A
i
s
~
V
n
i

t
~
V
n1
i
~
A
i
s
@
s
~
V
n
i


X
a
ik
2F
int
Z
a
ik
t
~
V
n1
i
~
P
k
ik
~
V
n
k

X
a
ik
2F
slip
Z
a
ik
t
~
V
n1
i
~
P
i
ik
H
ik
~
V
n
i
:
Then we deduce that X
n
i
+Y
n
i
= T
1
+T
2
+T
3
with
T
1
=
Z
T
i
X
s2fx;y;zg
t
@
s

~
Z
n1
i

p
A
i
0
1
~
A
i
s
p
A
i
0
1
~
Z
n
i

t
~
Z
n1
i
p
A
i
0
1
~
A
i
s
p
A
i
0
1
@
s

~
Z
n
i

;
T
2
= 
X
a
ik
2F
int
\@T
i
Z
a
ik
t
~
Z
n1
i
p
A
i
0
1
p
A
k
0
p
A
k
0
1
~
P
k
ik
p
A
k
0
1
~
Z
n
k
;
T
3
= 
X
a
ik
2F
slip
\@T
i
Z
a
ik
t
~
Z
n1
i
p
A
i
0
1
~
P
i
ik
p
A
i
0
1
p
A
i
0
H
ik
p
A
i
0
1
~
Z
n
i
:
The matrix
p
A
i
0
1
~
A
i
s
p
A
i
0
1
is symmetric and 

p
A
i
0
1
~
A
i
s
p
A
i
0
1

= 

A
i
s

,then we have
jT
1
j 
2
i

i
P
i
V
i
k
~
Z
n1
i
k
T
i
k
~
Z
n
i
k
T
i
=
2
i

i
P
i
V
i
kj
~
W
n1
i
kj
T
i
kj
~
W
n
i
kj
T
i
:
The matrix
p
A
k
0
1
~
P
k
ik
p
A
k
0
1
is also symmetric and 

p
A
k
0
1
~
P
k
ik
p
A
k
0
1

= 

P
k
ik

,then
jT
2
j 
X
a
ik
2F
int
\@T
i


P
k
ik

r


p
A
k
0
A
i
0
1
p
A
k
0

k
~
Z
n1
i
k
a
ik
k
~
Z
n
i
k
a
ik

X
a
ik
2F
int
\@T
i


P
k
ik

r


A
k
0
A
i
0
1

kj
~
W
n1
i
kj
a
ik
kj
~
W
n
i
kj
a
ik
:
Since the matrices A
i
0
and H
ik
commute and H
2
ik
= Id,we have
jT
3
j 
X
a
ik
2F
ref
\@T
i


P
i
ik

k
~
Z
n1
i
k
a
ik
k
~
Z
n
i
k
a
ik
=
X
a
ik
2F
ref
\@T
i


P
i
ik

kj
~
W
n1
i
kj
a
ik
kj
~
W
n
i
kj
a
ik
:
Using the above denitions and inequalities of the form ab 
1
2

a
2
+b
2

,rewriting F
n
i

1
2

kj
~
W
n
i
kj
2
T
i
+kj
~
W
n1
i
kj
2
T
i


t
2
(jT
1
j +T
2
j +T
3
j) leads to the intermediate result (5.14).
References
1.S.K.Lele,"Computational Aeroacoustics:a review",35th Aerospace Sciences Meeting & Exhibit,
Reno,January 6-10,AIAA Paper 97-0018 (1997).
A dissipation-free DGTD method for the three-dimensional linearized Euler equations
2.M.Lesieur and O.Metais,Annu.Rev.Fluid Mech.,28,pp.45-82 (1996).
3.C.K.W.Tam and J.C.Webb,J.Comput.Phys.,107,pp.262{281 (1993).
4.C.Bogey,C.Bailly and D.Juve,AIAA J.,40(2),pp.235-243 (2002).
5.C.K.W.Tam,AIAA J.,33(10),pp.1788-1796 (1995).
6.M.J.Lighthill,Proc.Roy.Soc.London,211,Ser.A,1107,pp.564-587 (1952)
7.G.M.Lilley,The generation and radiation of supersonic jet noise,Air Force Aero Propulsion
Laboratory,AFAPL-TR-72-53,Vol 4,Part 2 (1972).
8.C.Bogey,Calcul direct du bruit aerodynamique et validation de modeles acoustiques hybrides,
Ph.D.thesis,Ecole centrale de Lyon,2000.
9.L.Fezoui,S.Lanteri,S.Lohrengel and S.Piperno,RAIRO M2AN,39,(6)pp.1149-1176 (2005)
10.F.Q.Hu,J.Comput.Phys.,129,pp.201-219 (1996)
11.F.Q.Hu,J.Comput.Phys.,173,pp.455-480 (2001)
12.F.Q.Hu,M.Y.Hussaini and P.Rasetarinera,J.Comput.Phys.,151,pp.921-946 (1999)
13.D.Kroner,Math.Comput.,195(57) (1991)
14.E.Becache,S.D.Fauqueux and P.Joly INRIA Research Report RR-4304 (INRIA,2001).
15.C.K.W.Tam and Z.Dong,J.Comput.Phys.,4(2),pp.175-201 (1996).
16.S.Piperno and L.Fezoui,INRIA Research Report RR-4733 (INRIA,2003)
17.A.Harten,ICASE Research report,(81-34) (1981).
18.J.Hesthaven and T.Warburton,J.Comput.Phys.181 (1) 186{221 (2002).
19.P.Delorme and C.Peyret,Galerkin discontinous method for computational aeroacoustics,Private
communication.
20.M.Bernacki and S.Piperno,INRIA Research Report RR-4932 (INRIA,2003).
21.M.Bernacki and S.Piperno,INRIA Research Report RR-4699 (INRIA,2003).
22.M.Bernacki,S.Lanteri,and S.Piperno,to appear in J.Comput.Acoust..
23.C.K.W Tam and J.C.Hardin,Second computational aeroacoustics workshop on benchmark
problems,NASA Ohio Aerospace Institute (1997).
24.C.K.W Tam,J.C.Hardin and J.R.Ristorcelli,Workshop on benchmark problems in computa-
tional aeroacoustics,NASA Ohio aerospace institute (1995).
25.Fourth computational aeroacoustics workshop on benchmark problems NASA,Ohio Aerospace
Institute (2003).