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 satises a balance equation with

source term at the continuous level,and that our numerical framework satises 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

conrm 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.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

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 diculty of these approaches is the

dierence 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,ecient,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 veries 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 veries 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 veried.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 themarticially (like methods based on upwind

uxes,which damp instabilities,in an uncontrolled way though);

5.we show nally that,by introducing an articial 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 dierent congurations,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

Arst step towards aeroacoustic wave propagation consists in linearizing Euler equations

around any uniform ow dened by some constant physical elds (

0

;~v

0

;p

0

),therefore being

a steady-state solution of Euler equations.Innitely 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 dene a volumic

aeroacoustic energy E =

1

2

k

~

Wk

2

(a very simple expression in function of

~

W),which veries

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 dened 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 dicult 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 denite 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 dened by E =

1

2

t

~

WA

0

0

1

~

W

1

2

t

~

VA

0

0

~

V veries

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 articial 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 denitions 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 dene 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 dened 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 dene 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

= nt.Assuming we dispose of

~

W

n1

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

n1

i

2t

=

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 denite 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 veries 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 dicult 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

n1

i

+

~

W

n+1

i

2

,where

p

A

i

0

is the positive

square root of the symmetric denite 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

n1

i

.

Bernacki,Piperno

3.3.Energy balance and stability

In order to investigate stability,we dene 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

n1

i

+

t

4

X

a

ik

2F

abs

Z

a

ik

t

A

i

0

1

~

W

n1

i

M

ik

A

i

0

1

~

W

n1

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

n1

i

+

~

V

n+1

i

M

ik

~

V

n1

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 denite

formof numerical unknowns (

~

W

n1

i

;

~

W

n

i

) under some CFL-like sucient 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 coecients 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-diusive) 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 congurations,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 eciency.

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.Articial 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 eciency 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 eciency.

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,articial

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 eciency

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 eciency.

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 amplied 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 dierent 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.Dierent 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 prole

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 prole (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 prole,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 diculties 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 prole 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 prole,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 congurations.

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 satises a balance equa-

tion with source term at the continuous level,and that our numerical method satises

an equivalent balance equation at the discrete level.This shows that the method is gen-

uinely dissipation-free,which is conrmed by numerical results,including test-cases where

Kelvin-Helmholtz instabilities appear.

Further works can concern many dierent 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

n1

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

n1

i

M

ik

A

i

0

1

~

W

n+1

i

+

~

W

n1

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

n1

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

n1

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

n1

i

M

ik

~

V

n1

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

n1

i

+

~

V

n+1

i

M

ik

~

V

n1

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

n1

i

+

~

V

n+1

i

M

ik

~

V

n1

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 denite quadratic form of unknowns

~

W

n1

i

and

~

W

n

i

under the condition given in Eq.(3.12).Let us introduce the following denitions.

Denition 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 dene k

~

Xk

a

ik

= (

R

a

ik

k

~

Xk

2

)

1=2

.Finally,the matrix A

i

0

being positive denite,we dene a second norm by kj

~

Xkj

T

i

= k

p

A

i

0

1

~

Xk

T

i

.

Denition 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)

Denition 3 We choose here some notations.We dene

~

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 dene

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 denite quadratic form of unknowns

~

W

n1

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

n1

i

kj

2

T

i

i

P

i

i

t

V

i

kj

~

W

n

i

kj

T

i

kj

~

W

n1

i

kj

T

i

t

4

X

k2V

i

ik

ik

k~n

ik

k

V

i

kj

~

W

n1

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 denitions 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

n1

i

kj

2

T

i

i

i

t

V

i

kj

~

W

n

i

kj

T

i

kj

~

W

n1

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

n1

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

n1

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

n1

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

n1

i

kj

2

T

i

:

Therefore,the energy F

n

is a positive denite 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

n1

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

n1

i

M

ik

A

i

0

1

~

W

n1

i

+

~

W

n+1

i

=

1

2

X

a

ik

2F

abs

\@T

i

Z

a

ik

t

~

V

n1

i

M

ik

~

V

n1

i

+

~

V

n+1

i

Y

n

i

=

1

t

Z

T

i

t

~

W

n1

i

A

i

0

1

~

W

n+1

i

~

W

n1

i

= 2

Z

T

i

X

s2fx;y;zg

t

@

s

~

V

n1

i

~

A

i

s

~

V

n

i

X

k2V

i

Z

a

ik

t

~

V

n1

i

~

P

i

ik

~

V

n

i

+

~

P

k

ik

~

V

n

k

=

Z

T

i

X

s2fx;y;zg

t

@

s

~

V

n1

i

~

A

i

s

~

V

n

i

t

~

V

n1

i

~

A

i

s

@

s

~

V

n

i

X

k2V

i

Z

a

ik

t

~

V

n1

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

n1

i

~

A

i

s

@

s

~

V

n

i

X

a

ik

2F

int

Z

a

ik

t

~

V

n1

i

~

P

k

ik

~

V

n

k

X

a

ik

2F

slip

Z

a

ik

t

~

V

n1

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

n1

i

p

A

i

0

1

~

A

i

s

p

A

i

0

1

~

Z

n

i

t

~

Z

n1

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

n1

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

n1

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

n1

i

k

T

i

k

~

Z

n

i

k

T

i

=

2

i

i

P

i

V

i

kj

~

W

n1

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

n1

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

n1

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

n1

i

k

a

ik

k

~

Z

n

i

k

a

ik

=

X

a

ik

2F

ref

\@T

i

P

i

ik

kj

~

W

n1

i

kj

a

ik

kj

~

W

n

i

kj

a

ik

:

Using the above denitions 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

n1

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.Juve,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 aerodynamique et validation de modeles 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.Kroner,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).

## Σχόλια 0

Συνδεθείτε για να κοινοποιήσετε σχόλιο