Journal of Computational Acoustics, cfIMACS

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, [email protected] SERGE PIPERNO Project-team Caiman, CERMICS and INRIA, BP93, F06902 Sophia Antipolis Cedex, [email protected]

Received Revised

(to be inserted 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 solution. In the general context of a non-uniform supporting flow, we prove, using the well-known symmetrization of Euler equations, that some aeroacoustic energy satisfies a balance equation with source term at the continuous level, and that our numerical framework satisfies an equivalent balance equation at the discrete level and is genuinely dissipation-free. In the case of P1 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 confirm the theoretical properties of the method. They include test-cases where Kelvin-Helmholtz instabilities appear. Keywords: aeroacoustics; acoustic energy; linearized Euler equations; non-uniform steady-state flow; Discontinuous Galerkin method; time domain; energy-conservation.

1. Introduction Aeroacoustics is a domain where numerical simulation meets great expansion. The minimization of acoustic pollutions by aircrafts at landing and take off, or more generally by aerospace and terrestrial vehicles, is now an industrial concern, related to more and more severe norms. Different approaches coexist under the Computational Aeroacoustics activity. The most widely used methods belong to classical Computational Fluid Dynamics and consist in solving partial differential equations for the fluid, without distinction between the supporting (possibly steady-state) flow and acoustic perturbations 1 . The equations modeling the fluid can be Euler or Navier-Stokes equations, possibly including extended models like turbulence, LES techniques, etc 2 . One particular difficulty of these approaches is the difference in magnitude between the flow 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 determination of the supporting steady-state flow and in modeling the generation of noise (for example by providing equivalent acoustic sources), from the propagation of acoustic perturbations3,4,5 . For this problem, linearized Euler equations around the supporting flow are to be solved and provide a good description of the propagation of aeroacoustic perturbations 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 Lilley7 (a clear description can be found in a more recent reference 8 ). The noise source modeled or derived from the steady-state flow 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, efficient, 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 flows, obtained using a given Euler solver. The supporting flow 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 domain 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 equations16 . We use an element-centered formulation with centered numerical fluxes and an explicit leap-frog time scheme. This kind of method provides a dissipation-free approximation of propagation equations and allows for the accurate estimation of aeroacoustic energy variation, which is not possible with numerical methods (finite volumes, discontinuous Galerkin, spectral elements) based on upwind numerical fluxes. The fact that centered numerical fluxes in discontinuous Galerkin time-domain methods can lead to discretization methods inducing no numerical dissipation is quite well known. This was for example first 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 uniform supporting flow, at the continuous level (i.e. before space discretization), some quadratic energy verifies a balance equation without source term. This means energy is conserved (up to boundaries); 2. in this “uniform supporting flow” case, we are able to prove that our DG method (with leap-frog time-scheme and centered fluxes) 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 claim that we “control energy variations in the uniform case”; 3. accordingly, for a non-uniform supporting flow, at the continuous level (i.e. before space discretization), we use the well-known symmetrization of nonlinear Euler equations 17 to derive an aeroacoustic energy which verifies some balance equation with source term. Because of this unsigned source term, aeroacoustic waves can be damped or excited by the supporting flow. 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 uniform supporting flow” 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 term is also verified. 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 them artificially (like methods based on upwind fluxes, which damp instabilities, in an uncontrolled way though); 5. we show finally that, by introducing an artificial 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 uniform or non-uniform supporting flow, 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 different configurations, obtained with MPI-based Fortran parallel implementations 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 flow. Therefore, we linearize the three-dimensional Euler equations around a given steady flow and only take into account first-order perturbation terms. For a perfect inviscid gas, Euler equations read:

∂t

ρ ρu ρv ρw e

+ ∂x

ρu ρu2 + p ρuv ρuw (e + p)u

+ ∂y

ρv ρuv ρv 2 + p ρvw (e + p)v

+ ∂z

ρw ρuw ρvw ρw2 + p (e + p)w

= 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 − 12 ρk~v k2 ), where γ is a fixed constant (γ > 1). 2.1. Linearization around a uniform flow A first step towards aeroacoustic wave propagation consists in linearizing Euler equations around any uniform flow defined by some constant physical fields (ρ 0 , ~v0 , p0 ), therefore being a steady-state solution of Euler equations. Infinitely small perturbations (δρ, δ~v , δp) of this flow verify the following linear hyperbolic system of conservation laws: ~ + A x ∂x W ~ + A y ∂y W ~ + A z ∂z W ~ = ~0, ∂t W

(2.2)

t

~ = (δp − c2 δρ, t ρ0 c0 δ~v , δp), c0 corresponds to the uniform sound speed given by where W 0 2 ~ has been chosen such that the constant matrices A x , Ay , and c0 = γp0 /ρ0 . The variable W Az are symmetric. They are given by: u0 0 0 0 0 v0 0 0 0 0 w0 0 0 0 0 0 u0 0 0 c 0 0 v0 0 0 0 0 w0 0 0 0 Ax = 0 0 u0 0 0 , Ay=0 0 v0 0 c0 , Az= 0 0 w0 0 0 0 0 0 u0 0 0 0 0 v0 0 0 0 0 w 0 c0 0 c0 0 0 u 0 0 0 c 0 0 v0 0 0 0 c 0 w0 In the particular case of the linearization around a uniform flow, one can define a volumic ~ 2 (a very simple expression in function of W), ~ which verifies aeroacoustic energy E = 21 kWk ~ = 0, where the energy flux F ~ is given by the balance equation ∂t E + divF

1t ~ ~ s ∈ {x, y, z}. WAs W, 2 Thus, away from boundaries, the aeroacoustic energy E is exactly conserved. Fs =

2.2. Linearization around a non-uniform flow A more interesting framework consists in linearizing Euler equations around a nonuniform steady-state solution. In that case, the steady flow is defined by smoothly varying physical quantities (ρ0 , ~v0 , p0 ). Linearizing straightforwardly Euler equations (2.1) yields: ~ + ∂x A0x W ~ + ∂y A0y W ~ + ∂z A0z W ~ = 0, ∂t W (2.3)

~ now denotes the perturbations of conservative variables (i.e. W ~ T = (δρ, ρ0 δ~v + where W ~v0 δρ, δp/(γ − 1) + ρ0~v0 .δ~v + k~v0 k2 δρ/2)) and the space-varying matrices A 0x , A0y , and A0z are given in function of γ˜ = γ − 1, α0 = c20 /˜ γ + k~v0 k2 /2, β0 = (γ − 2)kV0 k2 /2 − c20 /˜ γ , and 3 the canonical basis (~ex , ~ey , ~ez ) of R by t 0 ~es 0 A0s = γ2˜ k~v0 k2~es − (~v0 .~es )~v0 (~v0 .~es )I3 − γ˜~es t~v0 + ~v0 t~es γ˜~es , s ∈ {x, y, z}. (2.4) t t β0 (~v0 .~es ) α0 ~es − γ˜ (~v0 .~es ) ~v0 γ(~v0 .~es )

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

In this equation, the matrices A0x , A0y , and A0z are not symmetric anymore, and it is very difficult to deduce any aeroacoustic energy balance equation. Therefore, we consider other acoustic variables, derived from the quite classical symmetrization of Euler equations. Assuming the flow is smooth enough, the change of variables (ρ, ρ~v , e) → γ γ e˜ γ v , − ρ˜ (− p + γ + 1 − ln ρpγ , ρ˜ p~ p ) transforms Euler equations (2.1) into a symmetric system of conservation laws (i.e. Jacobians of fluxes 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 ~ + ∂x A ˜ 0V ~ + ∂y A ˜ 0V ~ + ∂z A ˜ 0V ~ = 0, A00 ∂t V (2.5) x y z ~ is given in function of variables W ~ by V ~ = A0 −1 W ~ and where V 0 t 1 ~v0 α0 − c20 /γ ρ0 ˜0 c20 A00 = ~v0 v0 t~v0 α0~v0 ; As = A0s A00 , s ∈ {x, y, z}. (2.6) γ I3 + ~ γ˜ α0 − c20 /γ α0 t~v0 α20 − c40 / (γ˜ γ)

A00 clearly is symmetric and it can be proved that it is definite positive (and then not ~ by A0 V ~ in Eq. 2.3 (and by singular). Eq. 2.5 can also be obtained simply by replacing W 0 0 ˜s noting that ∂t A0 = 0). Finally, the reader can also check that the symmetric matrices A (s ∈ {x, y, z}) are given by t 0 ~es (~v0 .~es ) ˜ 0 = (~v0 .~es ) A0 + p0 ~es (~v0 .~es )~v0 + α0~es . ~es t~v0 + ~v0 t~es A s 0 γ˜ t t (~v0 .~es ) (~v0 .~es ) ~v0 + α0 ~es 2α0 (~v0 .~es )

t ~ 0 −1 W ~ ≡ 1 t VA ~ 0V ~ verifies Then, the volumic aeroacoustic energy E defined by E = 21 WA 0 0 2 the following balance equation with source term: ( t ~A ˜ 0 V, ~ s ∈ {x, y, z}. Fs = V sh i ~ ∂t E + divF = S, with (2.7) t ~ ~ ∂x (A ˜ 0 ) + ∂ y (A ˜ 0 ) + ∂ z (A ˜ 0 ) V. S = − 21 V x y z

Thus the aeroacoustic energy is not conserved and the variations in the steady flow considered can damp or amplify aeroacoustic waves, unless the source term vanishes (which is the case for a uniform flow for example). In the sequel, we shall mainly discretize the conservative 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 flow 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 field artificial boundary ∂Ω∞ ). Ω is partitioned into a finite number of polyhedra (each one having a finite number of faces). For each polyhedron T i , called ”control volume” or ”cell”, Vi 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 T T F int = F/∂Ω, physical faces F phys = F ∂Ωphys and absorbing faces F abs = F ∂Ω∞ . For T each internal face aik = Ti Tk , we denote by Sik the measure of aik and by ~nik the unitary normal, oriented from Ti towards Tk . The same definitions are extended to boundary faces, the index k corresponding to a fictitious 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 P face in common). We also define the perimeter P i of Ti by Pi = k∈Vi Sik . We recall the P nik = 0. following geometrical property for all control volumes: k∈Vi Sik ~ Following the general principle of discontinuous Galerkin finite element methods, the unknown field inside each control volume is seeked for as a linear combination of local basis vector fields ϕ ~ ij , 1 ≤ j ≤ di (generating the local space Pi ) and the approximate field is allowed to be fully discontinuous across element boundaries. Thus, a numerical flux function has to be defined to approximate fluxes at control volumes interfaces, where the approximate solution is discontinuous. This context is quite general. Actual implementations of the method have been considered only on tetrahedral meshes, where control volumes are the tetrahedra themselves. We shall only consider constant (P0 ) or linear (P1 ) approximations inside tetrahedra. 3.1. Time and space discretizations We only consider here the most general case of aeroacoustic wave propagation in a nonuniform steady flow. Also, in order to limit the amount of computations, we restrict our study to piecewise constant matrices A 0s (s ∈ {x, y, z}) given in Eq. (2.4). For each control volume Ti , for s ∈ {x, y, z}, we denote by Ais an approximate for the average value of A 0s over Ti . Dot-multiplying Eq. (2.3) by any given vector field ϕ ~ , integrating over T i and integrating by parts yields Z

Z Z X X ~ ∂W t ~ − ~ . ϕ ~· ∂s ϕ ~ A0s W = ϕ ~· ns A0s W ∂t Ti Ti ∂Ti s∈{x,y,z}

(3.8)

s∈{x,y,z}

~ by the approximate field W ~ i and Inside volume integrals over Ti , we replace the field W ~ the matrices A0s by their respective average values A is . For boundary integrals over ∂Ti , W

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

is discontinuous, and we define totally centered numerical fluxes, i.e.: h i ( ~ i + Pk W ~ k , ~ ∀i, ∀k ∈ Vi , (nikx A0x + niky A0y + nikz A0z )W ' 21 Piik W ik |aik

with Piik = nikx Aix + niky Aiy + nikz Aiz ,

Pkik = nikx Akx + niky Aky + nikz Akz .

(3.9)

Concerning the time discretization, we use a three-level leap-frog scheme. The unknowns ~ i are approximated at integer time-stations t n = n∆t. Assuming we dispose of W ~ n−1 W i n+1 ~ n , the unknowns W ~ and W are seeked for in Pi such that, ∀~ ϕ ∈ Pi , i i Z

Ti

ϕ ~·

~ n+1 − W ~ n−1 Z W i i = 2∆t Ti

X

t

~ in − ∂s ϕ ~ Ais W

s∈{x,y,z}

XZ

ϕ ~·

aik

k∈Vi

~ n + Pk W ~ n Piik W i ik k . 2

(3.10)

Again, the time scheme above is almost explicit. Each time step only requires the inversion of local symmetric positive definite matrices of size (d i × di ). In the particular case where Pi is a complete linear (P1 ) 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 flow and the acoustic perturbations. This means that we assume that for any slip boundary face a ik belonging to the control volume Ti , the steady solution of Euler equations verifies a slip condition at the discrete level, i.e. ~nik · ~vi0 = 0. For the acoustic perturbations, we use a mirror fic~ k in the computation of the boundary flux given in Eq. (3.10). We take titious state W δρk = δρi , δpk = δpi , and δ~vk = δ~vi − 2(~nik · ~vi )~nik (which implies (δ~vk − δ~vi ) × ~nik = 0 and δ~vk .~nik = −δ~vi .~nik ). For an absorbing boundary face aik , upwinding is used to select outgoing waves only. Before discretization in time, classical upwinding leads to a boundary flux F ik given by + i ~ Fik = Pik Wi , where for any diagonalizable matrix Q = S −1 DS with D diagonal, Q+ = (Q + |Q|)/2 and terms of the diagonal matrix |D| are the moduli of the eigenvalues. This ~ i . However, for this intuitive numerical flux, it is ~ k = |Pi |W general idea leads to Pkik W ik very difficult to prove that the resulting time-scheme is stable and that energy is actually sent in the exterior domain. We then consider the numerical flux based on the following √ √ i √ i −1 i √ i √ i −1 W ~ n−1 +W ~ n+1 n k i i ~ , where Ai is the positive fictitious state: P W = A A P A A ik

k

0

0

ik

0

0

2

0

square root of the symmetric definite positive matrix A i0 . Indeed, this expression derives from the intuitive upwind flux 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 √ i √ i −1 i √ i √ i −1 n−1 n 22 k ~ = A A P A A W ~ available . It takes the form Pik W . 0 0 0 0 i k ik

Bernacki, Piperno

3.3. Energy balance and stability In order to investigate stability, we define a discrete aeroacoustic energy F n by: Z t 1X n ~ n Ai −1 W ~ n + tW ~ n+1 Ai −1 W ~ n−1 W F = i 0 i 0 i i 2 T i i Z t ∆t X −1 ~ n−1 i −1 ~ n−1 n+1 ~ + , M A W + W Ai0 W ik 0 i i i 4 abs aik aik ∈F

√ √ −1 i √ i −1 √ i √ √ −1 √ √ ˜ A A . One can show that the with Mik = Ai0 Ai0 Piik Ai0 Ai0 ≡ Ai0 Ai0 P 0 0 ik matrices Mik 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 − Fn ): Z t n k ∆t X n+1 n ~ (P ˜ −P ˜ i )V ~ n+1 + t V ~ n+1 (P ˜k − P ˜ i )V ~n V F −F = − i ik ik ik ik k i k 2 aik aik ∈F int Z t ∆t X ~ n−1 + V ~ n+1 Mik V ~ n−1 + V ~ n+1 . − V (3.11) i i i i 4 abs aik aik ∈F

~ n ≡ Ai −1 W ~ n , ∀i, ∀n. The first term is a where we have used the auxiliary variables V 0 i i 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 influence on the global energy balance. In order to prove stability, we show in Annex 2 that F n is a quadratic positive definite ~ n ) under some CFL-like sufficient stability condition ~ n−1 , W form of numerical unknowns (W i i on the time-step ∆t: 2Vi ∀i, ∀k ∈ Vi , ∆t (2λi αi + βik ρik ) < , (3.12) Pi where αi and βik are dimensionless regularity coefficients depending of basis functions and element aspectratio, λi = |ui0 | + |v0i| + |w0i | +3ci0 , and ρik = |~v0i ·~nik | + ci0 fora boundary face 2 2 −1 −1 for an internal , |~v0k · ~nik | + ck0 ρ Ai0 Ak0 and ρ2ik = sup |~v0i · ~nik | + ci0 ) ρ Ak0 Ai0 face (ρ here denotes the spectral radius of a matrix). In the case of a uniform flow, we have P iik = Pkik . Thus the aeroacoustic energy is non-increasing (and exactly conserved if no absorbing boundary is present, which shows the scheme is genuinely non-diffusive) and the scheme is stable under a CFL-type stability condition depending on the size of elements and sup i (k~vi k + c0i ). 4. Numerical results We dispose of a three-dimensional parallel implementation of the DGTD method presented in the previous section. Any subsonic steady flow can be considered, even with strong flow gradients. However, the flow, given as the output of a non-linear Euler equations solver, has to be post-processed: average of the flow 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 testcases in two and three space dimensions, in order to validate the method on benchmark problems, test the method on complex flows and configurations, and finally evaluate the performance of the parallel Fortran 77 implementation, based on the MPICH implementation 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 efficiency. 4.1. Acoustic perturbations in a three-dimensional uniform flow This quite classical test-case23,24 is aimed at validating the DGTD method proposed by mixing acoustic, entropic and vorticity perturbations, in a uniform horizontal flow (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 M1 M2

# vertices 68,921 531,441

# tetrahedra 384,000 3,072,000

# absorbing faces 19,200 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 finer mesh M2. Artificial reflections 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 flow is uniform, the total acoustic energy is first exactly conserved, and then decreasing, once perturbations have reached the absorbing boundaries (at t = 15). Computation times are given in Table 2. Parallel efficiency reaches a satisfying level and total CPU time are really acceptable. 4.2. Acoustic perturbations in a three-dimensional subsonic Couette flow This second test-case is more complex and is classically used to evaluate the merits of a numerical method for non-uniform shear flows 25 . The steady flow considered is horizontal, with ~v0 = 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 flow. 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.

3.5 3

F

2.5 2

1.5 1 0.5

0

5

10 15 20 25 30 35 40 45 50 55 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 efficiency. Mesh M1 M2 -

Np 1 4 8 16 8 16

CPU 33381 s 8167 s 4286 s 2195 s 34012 s 17203 s

REAL 33447 s 8412 s 4465 s 2381 s 35651 s 18760 s

% CPU 100% 97% 96% 92% 95% 92%

S(Np ) 1 4 7.5 15.2 1 1.98

waveguide, the aim of this test-case is to validate the accurate propagation and convection of the waves by the shear flow, 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 flow test-case. Characteristics of meshes. Mesh M1 M2

# vertices 73,629 522,801

# tetrahedra 405,600 3,000,000

# absorbing faces 23,504 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 flow : in the numerical results, artificial reflections from absorbing boundaries have a relative amplitude of 1.5% for the coarse mesh M1 (respectively 1% for the finer 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 flow (thus the bottom part of this perturbation travels slowly with the supporting flow, whose velocity vanishes in the plane z = 0). The time evolution of the global acoustic energy is given on Fig. 4. As expected, the energy is not conserved anymore because the non-uniform flow provides energy to the acoustic perturbations. The energy decreases when perturbations reach the absorbing boundaries. A part of the energy corresponding to the bottom part of the entropic perturbation remains in the domain. Computation times are given in Table 4 and reveal similar parallel efficiency 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 flows. 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 flow 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.08 0.07 0.06

Energy

0.05 0.04 0.03 0.02 0.01 0

0

1

2

3 Time (s)

4

5

Fig. 4. Couette flow test-case. Time evolution of the acoustic energy.

Table 4. Couette flow test-case. Parallel CPU times and efficiency. Mesh M1 M2 -

Np 4 8 16 8 16

CPU 12830 s 6428 s 3230 s 46359 s 23412

REAL 13008 s 6590 s 3373 s 47743 s 24959 s

% CPU 99% 98% 96% 97% 94%

S(Np ) 1 2 3.9 1 1.98

Bernacki, Piperno

acoustic energy in the flow. 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 flow 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, p0 = 1/γ. In this flow, u0 goes through an inflection point at y = 0.25 (and k~v 0 k = 0.5 on the inflection point). In this particular flow, the acoustic energy is amplified along the inflection line and at the same time, perturbations are convected along this line, which leads to an instability. The behavior is completely different from the one obtained with a shear flow of the form ~v0 = (y + 0.25)~ex , 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 flow with an inflection point. Different means are under current investigation to control this kind of instability.

Fig. 5. Contours of kδ~v k at t=0.257 for the linear (left) and inflection (right) shear flows with source term

4.4. Two-dimensional instabilities past a NACA airfoil profile This test-case is aimed at evaluating the occurrence of Kevin-Helmholtz-type instabilities for steady-state flows with sharp gradients. We have considered a test-case past a NACAtype airfoil profile (courtesy of ONERA 19 ). The steady-state solution at M∞ = 0.5 obtained by our finite-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 profile, and our absorbing boundary conditions in the far field boundary. For this testcase, 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δ~v k are shown on Fig. 7. Some numerical treatments are possible to overcome these difficulties and get rid of instabilities. One of them4 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δ~v k are shown on the right part of Fig. 8 and can be compared to classical acoustic scattering (uniform still flow). The comparison of contours near the trailing edge for acoustics and aeroacoustics is shown on Fig. 9, where the shear flow at the trailing edge has a visible influence 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 flow past the airfoil profile at M∞ = 0.5

Fig. 7. Contours of kδ~v k at t = 1.57ms for the original aeroacoustic model.

Bernacki, Piperno

propagation of acoustics (and generation of vortices).

Fig. 8. Contours of kδ~v k for the stabilized aeroacoustic model (left) and classical acoustics (right).

Fig. 9. Tailing edge zoom – kδ~v k 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 flow past a sphere and the propagation of waves emitted by an acoustic source (Gaussian profile, periodic source term with period T = 0.2s. The subsonic flow with M ∞ = 0.5 is obtained by a three-dimensional MUSCL-extended finite-volume solver of three-dimensional Euler equations on tetrahedral grids. The steady-state solution presents Mach numbers between 0.002 and 0.8 and contours 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 field 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 324,471

# elements 1,870,288

# abs. faces 10,092

# slip faces 46,080

Np 16

CPU 50h4mn

REAL 51h18mn

% CPU 97,6%

resulting propagated waves for this aeroacoustic case in a steady flow past a sphere (after 10 and 25 periods). The numerical results (scattering and reflection 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 configurations. 5. Conclusion In this paper, we have presented a time domain Discontinuous Galerkin dissipationfree 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 flow, we have proved, using the well-known symmetrization of Euler equations, that the aeroacoustic energy satisfies a balance equation with source term at the continuous level, and that our numerical method satisfies an equivalent balance equation at the discrete level. This shows that the method is genuinely dissipation-free, which is confirmed by numerical results, including test-cases where Kelvin-Helmholtz instabilities appear. Further works can concern many different 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 flows 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 flow (currently, it is only P0 ). 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 i0 is symmetric, we have

F

n+1

−F

n

=

1X 2 i

∆t + 4

Z

t

Ti

~ n+1 Ai −1 W ~ n+1 − W ~ n−1 ~ n Ai −1 W ~ n+2 − W ~ n + tW W 0 i 0 i i i i i

X

aik ∈F abs

Z

aik

t

~ n+2 + W ~ n ~ n Mik Ai −1 W W 0 i i i

Ai0

−1

t

−1 ~ n−1 Ai0 W i

−

−1 Mik Ai0

~ n+1 + W ~ n−1 W i i

.

The variational form of the scheme (3.10) allows to consider any field in P i , thus for example

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

~ n−1 or W ~ n themselves. We obtain that Fn+1 = Fn + ∆t(T1 − 1 T2 − 1 T3 ) with the fields W i i 2 2 XZ

T1 =

i

Ti s∈{x,y,z}

X Z

T2 =

aik ∈F int

aik ∈F int

T3 =

X

aik ∈F slip

+

" X Z t

aik ∈F abs

aik

Z

t

aik

t

aik

X Z

+

t

X

t ˜iV ~n ~ n+1 A ˜iV ~ n+1 + ∂s V ~n A ∂s V s i s i i i

~ in P ˜i V ~ n+1 + P ˜k V ~ n+1 + t V ~ n+1 + P ~ n+1 ˜k V ˜i V ~n P V ik i ik k ki k ki i k t

aik

t ~ n+1 P ˜k V ~ n ˜i ~ n ~ n+1 P ˜i V ~ n ˜k ~ n V ki k + Pki Vi ik i + Pik Vk + Vk i

˜ i Hik V ~n ˜i V ~n+P ~ n+1 P ˜i V ~ n+1 + P ˜ i Hik V ~ n+1 + t V ~n P V ik i ik i ik i ik i i i

! !# ~ n+ V ~ n+2 ~ n−1 + V ~ n+1 V V t i i i ˜i V ˜i V ~ n+1 + Mik i ~ n+1 P ~n P + V ik i ik i + Mik i 2 2 Z t n 1 X ~ n−1 Mik V ~ n−1 + V ~ n+1 . ~ Mik V ~n+V ~ n+2 − t V − V i i i i i i 2 abs aik

~ in V

aik ∈F

˜ is symmetric, and since P ˜ i = −P ˜ i and P ˜ k = −P ˜ k , we obtain that Since P ik ki ik ki X Z

T2 = 2

aik

aik ∈F int

+

t

X Z

aik

aik ∈F int

T2 + T 3 = 2

XXZ i

+

k∈Vi

aik ∈F int

+

1 2

t aik

X Z X

t

aik ∈F abs

Z

~ n+1 , and ˜i V ˜k − P ~ n+1 + t V ~n P ˜i V ˜k − P ~n P V i ik ik k ik ik i k

~ nP ˜ i ~ n+1 + V i ik Vi

X

aik ∈F slip t

aik

~ nP ˜ i ~ n+1 + t V ~ nP ˜ k ~ n+1 V i ik Vi k ki Vk

Z

t

aik

~ nP ˜i ~ n+1 + t V ~ n+1 P ˜ i Hik V ~n V i ik Hik Vi ik i i

~ n+1 ~ n+1 + t V ˜i V ˜i V ~n P ˜k − P ~n P ˜k − P V i ik ik ik ik k i k t

aik

~ n−1 + V ~ n+1 Mik V ~ n−1 + V ~ n+1 . V i i i i

In the expressions above, the matrix H ik denotes the matrix used for slip boundary con~ n is given by W ~ n = Hik W ~ n (Hik does not ditions, in the sense that the fictitious state W i k k change the density and energy perturbations, and operates an orthogonal symmetry on the velocity perturbation). One can also show that H 2ik = I and that Hik commutes with Ai0 .

Bernacki, Piperno

We finally get that Fn+1 = Fn + ∆t(T4 + T5 + T6 ) with XZ X Z X t t t n i ~ P ˜ ~ n+1 ˜iV ~n − ˜iV ~ n+1 + ∂s V ~ n+1 A ~n A V T4 = ∂s V i ik Vi s i s i i i Ti s∈{x,y,z}

i

T5 = −

1 2

T6 = −

1 4

−

1 2

X Z

t

aik

aik ∈F ref

~ in P ˜ i Hik + V ik

Z

t

X Z

t

X

aik ∈F abs

aik ∈F int

aik

aik

k∈Vi

t

˜ i Hik P ik

aik

~ n+1 V i

~ n−1 + V ~ n+1 Mik V ~ n−1 + V ~ n+1 V i i i i

˜k − P ˜i V ˜k − P ˜i V ~ n+1 ~ n+1 + t V ~n P ~ in P V ik ik k ik ik i k

˜ s are symmetric and constant over control volumes, we get by integratSince the matrices A ˜ i Hik is skew-symmetric, we get T5 = 0. ing by parts that T4 = 0. Also, since the matrix P ik ~ variables, we get the result announced. Using the expression of Mik and using the W

Annex 2 ~ n−1 and W ~ n We shall prove that Fn is a positive definite quadratic form of unknowns W i i under the condition given in Eq. (3.12). Let us introduce the following definitions. R ~ ∈ Pi , we denote by kXk ~ T = ( kXk ~ 2 )1/2 the L2 -norm of X ~ over Ti . Definition 1 ∀X i Ti R 2 1/2 ~ ~ a =( . Finally, the matrix Ai0 For any interface aik ∈ F, we also define kXk ik aik kXk ) √ ~ T. ~ T = k Ai −1 Xk being positive definite, we define a second norm by k| Xk| i i 0 Definition 2 We assume some regularity for the basis functions ϕ ~ ij , 1 ≤ j ≤ di . More precisely, we assume that, for each control volume T i , there exist dimensionless positive constants αi and βik (for each k ∈ Vi ) such that ( αi Pi ~ ~ ~ ∈ Pi , ∀s ∈ {x, y, z}, k∂s XkTi ≤ Vi kXkTi , ∀X (5.13) ~ 2. ~ 2 ≤ βik k~nik k kXk ∀aik ∈ F, kXk aik

Vi

Ti

√ ~ n = Ai −1 W ~ n . We introduce Definition 3 We choose here some notations. We define Z 0 i i P i i i i i λ = s ρ As (and we recall ρ As = |~es .~v0 | + c0 ). For boundary faces aik , we recall that some neighboring fictitious control volume T k is imagined, and we take by convention: ~ n k|T , βki = βik , Vk = Vi . Finally, we define ρik for ~ n k|T = k|W ~ n kT , k|W ~ n kT = k W kW i i i k i k k k any face aik by: 2 2 −1 −1 , |~v0k · ~nik | + ck0 ρ Ai0 Ak0 aik ∈ F int , ρ2ik = sup |~v0i · ~nik | + ci0 ρ Ak0 Ai0 aik ∈ F slip , ρik = ρ Piik = |~v0i · ~nik | + ci0 aik ∈ F abs , ρik = 0,

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

~ n−1 and W ~ n In order to prove that Fn is a positive definite quadratic form of unknowns W i i under the condition (3.12), we first show the following result: i 1 ~ n 2 ~ n−1 k|2 − αi Pi λ ∆t k|W ~ n k|T k|W ~ n−1 k|T k|Wi k|Ti + k|W Fni ≥ Ti i i i i i 2 Vi ∆t X ρik βik k~nik k ~ n−1 2 ρki βki k~nki k ~ n 2 − (5.14) k|Wi k|Ti + k|Wk k|Tk . 4 Vi Vk k∈Vi

We can show easily that the above lemma gives the result. Using the above definitions and P ρki ∆t ~ n 2 assumptions, we have Fni ≥ k∈Vi k~nik k(Xik − βki4V k|Wk k|Tk ) with k i 1 βik ρik ∆t 1 ~ n 2 ~ n k|T k|W ~ n−1 k|T ~ n−1 k|2 − αi λ ∆t k|W k|Wi k|Ti + − k|W Xik = i Ti i i i i 2Pi 2Pi 4Vi Vi 1 αi λi ∆t 1 βik ρik ∆t αi λi ∆t ~ n k|2 + ~ n−1 k|2 − − − ≥ k|W k|W i Ti Ti i 2Pi 2Vi 2Pi 4Vi 2Vi X X k~nik kZik with k~nik kYik + Finally, Fn ≥ aik ∈F int

Yik

Zik

aik ∈F slip

S

F abs

1 βik ρik ∆t αi λi ∆t ~ n 2 ~ n−1 k|2 − − k|Wi k|Ti + k|W = Ti i 2Pi 4Vi 2Vi βki ρki ∆t αk λk ∆t ~ n 2 1 ~ n−1 k|2 , − − k|Wk k|Tk + k|W + Tk k 2Pk 4Vk 2Vk βik ρik ∆t αi λi ∆t ~ n 2 1 ~ n−1 k|2 . = − − k|Wi k|Ti + k|W Ti i 2Pi 4Vi 2Vi

Therefore, the energy Fn is a positive definite quadratic form of unknowns under the sufficient CFL-type stability condition proposed in Eq. (3.12). End ofthe proof. We now give the proof for the preliminary result (5.14). We first have n−1 2 1 n 2 n n n ~ ~ Fi = 2 k|Wi k|Ti + k|Wi k|Ti + ∆t 2 (Xi + Yi ) with Z X t 1 −1 ~ n−1 i −1 ~ n−1 n+1 n ~ M A W + W Ai0 W Xi = ik 0 i i i 2 aik aik ∈F abs ∩∂Ti Z X t n−1 1 ~ ~ n−1 + V ~ n+1 = V M V ik i i i 2 aik aik ∈F abs ∩∂Ti Z t 1 n ~ n−1 Ai −1 W ~ n+1 − W ~ n−1 W Yi = 0 i i i ∆t Ti Z X t XZ t n−1 i n ˜ V ~ − ~ ˜i V ~n+P ˜k V ~n ~ n−1 P = 2 A ∂s V V i

Ti s∈{x,y,z}

=

Z

X

Ti s∈{x,y,z}

s

i

k∈Vi

t

aik

i

ik

i

ik

XZ t n−1 i n n ˜ ∂s V ˜iV ~ ~ ~ ~ n−1 A − A − V ∂s V s i s i i i k∈Vi

aik

t

k

˜k V ~n ~ n−1 P V ik k i

Bernacki, Piperno

Using Hik and eliminating terms corresponding to absorbing boundaries leads to Z X t t n−1 i i n n n1 n n ˜ V ˜ ∂s V ~ ~ ~i ~i A A Xi + Y i = ∂s V s i − Vi s Ti s∈{x,y,z}

−

X Z

aik ∈F int

t aik

~ n−1 P ˜k V ~n V ik k − i

X

aik ∈F slip

Z

t aik

~ n. ~ n−1 P ˜ i Hik V V i ik i

Then we deduce that Xni + Yni = T1 + T2 + T3 with Z √ X t √ √ √ ~ n−1 Ai −1 A ˜ i Ai −1 Z ~ n − tZ ~ n−1 Ai −1 A ˜ i Ai −1 ∂s Z ~n , ∂s Z T1 = 0 s 0 i 0 s 0 i i i Ti s∈{x,y,z}

T2 = − T3 = −

Z

X

aik ∈F int ∩∂Ti

X

aik ∈F slip ∩∂Ti

t aik

Z

aik

√ √ √ √ ˜ k Ak −1 Z ~ n, ~ n−1 Ai −1 Ak Ak −1 P Z 0 0 0 ik 0 k i

t

√ √ √ √ ˜ i Ai −1 Ai Hik Ai −1 Z ~ ni . ~ n−1 Ai0 −1 P Z 0 0 0 ik i

√ −1 √ −1 √ −1 i √ i −1 ˜ A ˜ i Ai = ρ Ais , then we have The matrix Ai0 A is symmetric and ρ Ai0 A s s 0 0

2λi αi Pi ~ n−1 ~ n kT = 2λi αi Pi k|W ~ n−1 k|T k|W ~ n k|T . k Z i k Ti k Z i i i i i i Vi Vi √ −1 √ −1 √ −1 k √ k −1 ˜ A ˜ k Ak = ρ Pkik , then is also symmetric and ρ Ak0 P The matrix Ak0 P 0 0 ik ik |T1 | ≤

|T2 | ≤ ≤

ρ

Pkik

aik ∈F int ∩∂Ti

r √ −1 √ k ~ n ka ~ n−1 ka kZ ρ Ak0 Ai0 A0 kZ i ik ik i

X

ρ

Pkik

r −1 ~ in k|a . ~ n−1 k|a k|W k|W ρ Ak0 Ai0 ik ik i

X

aik

∈F int ∩∂T

i

Since the matrices Ai0 and Hik commute and H2ik = Id, we have X X n−1 ~ n k|a . ~ n−1 k|a k|W ~ n ka = ~ ρ Piik k|W kaik kZ |T3 | ≤ ρ Piik kZ i i ik ik ik i i aik ∈F ref ∩∂Ti

aik ∈F ref ∩∂Ti

1 2 + b2 , rewriting Fn ≥ Using the above definitions and inequalities of the form ab ≤ a i 2 n−1 2 ∆t 1 2 n ~ ~ k|Ti − 2 (|T1 | + T2 | + T3 |) leads to the intermediate result (5.14). 2 k|Wi k|Ti + k|Wi 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. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.

M. Lesieur and O. Metais, Annu. Rev. Fluid Mech., 28, pp.45-82 (1996). C.K.W. Tam and J.C. Webb, J. Comput. Phys., 107, pp.262–281 (1993). C. Bogey, C. Bailly and D. Juv´e, AIAA J., 40(2), pp.235-243 (2002). C.K.W. Tam, AIAA J., 33(10), pp.1788-1796 (1995). M.J. Lighthill, Proc.Roy. Soc. London, 211, Ser.A, 1107, pp.564-587 (1952) 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). C. Bogey, Calcul direct du bruit a´erodynamique et validation de mod`eles acoustiques hybrides, Ph.D. thesis, Ecole centrale de Lyon, 2000. L. Fezoui, S. Lanteri, S. Lohrengel and S. Piperno, RAIRO M2AN, 39, (6)pp.1149-1176 (2005) F.Q. Hu, J. Comput. Phys., 129, pp.201-219 (1996) F.Q. Hu, J. Comput. Phys., 173, pp.455-480 (2001) F.Q. Hu, M.Y. Hussaini and P. Rasetarinera, J. Comput. Phys., 151, pp.921-946 (1999) D. Kr¨ oner, Math. Comput., 195(57) (1991) E. Becache, S.D. Fauqueux and P. Joly INRIA Research Report RR-4304 (INRIA, 2001). C.K.W. Tam and Z. Dong, J. Comput. Phys., 4(2), pp.175-201 (1996). S. Piperno and L. Fezoui, INRIA Research Report RR-4733 (INRIA, 2003) A. Harten, ICASE Research report, (81-34) (1981). J. Hesthaven and T. Warburton, J. Comput. Phys. 181 (1) 186–221 (2002). P. Delorme and C. Peyret, Galerkin discontinous method for computational aeroacoustics, Private communication. M. Bernacki and S. Piperno, INRIA Research Report RR-4932 (INRIA, 2003). M. Bernacki and S. Piperno, INRIA Research Report RR-4699 (INRIA, 2003). M. Bernacki, S. Lanteri, and S. Piperno, to appear in J. Comput. Acoust.. C.K.W Tam and J.C. Hardin, Second computational aeroacoustics workshop on benchmark problems, NASA Ohio Aerospace Institute (1997). C.K.W Tam, J.C. Hardin and J.R. Ristorcelli, Workshop on benchmark problems in computational aeroacoustics, NASA Ohio aerospace institute (1995). Fourth computational aeroacoustics workshop on benchmark problems NASA, Ohio Aerospace Institute (2003).

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, [email protected] SERGE PIPERNO Project-team Caiman, CERMICS and INRIA, BP93, F06902 Sophia Antipolis Cedex, [email protected]

Received Revised

(to be inserted 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 solution. In the general context of a non-uniform supporting flow, we prove, using the well-known symmetrization of Euler equations, that some aeroacoustic energy satisfies a balance equation with source term at the continuous level, and that our numerical framework satisfies an equivalent balance equation at the discrete level and is genuinely dissipation-free. In the case of P1 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 confirm the theoretical properties of the method. They include test-cases where Kelvin-Helmholtz instabilities appear. Keywords: aeroacoustics; acoustic energy; linearized Euler equations; non-uniform steady-state flow; Discontinuous Galerkin method; time domain; energy-conservation.

1. Introduction Aeroacoustics is a domain where numerical simulation meets great expansion. The minimization of acoustic pollutions by aircrafts at landing and take off, or more generally by aerospace and terrestrial vehicles, is now an industrial concern, related to more and more severe norms. Different approaches coexist under the Computational Aeroacoustics activity. The most widely used methods belong to classical Computational Fluid Dynamics and consist in solving partial differential equations for the fluid, without distinction between the supporting (possibly steady-state) flow and acoustic perturbations 1 . The equations modeling the fluid can be Euler or Navier-Stokes equations, possibly including extended models like turbulence, LES techniques, etc 2 . One particular difficulty of these approaches is the difference in magnitude between the flow 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 determination of the supporting steady-state flow and in modeling the generation of noise (for example by providing equivalent acoustic sources), from the propagation of acoustic perturbations3,4,5 . For this problem, linearized Euler equations around the supporting flow are to be solved and provide a good description of the propagation of aeroacoustic perturbations 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 Lilley7 (a clear description can be found in a more recent reference 8 ). The noise source modeled or derived from the steady-state flow 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, efficient, 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 flows, obtained using a given Euler solver. The supporting flow 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 domain 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 equations16 . We use an element-centered formulation with centered numerical fluxes and an explicit leap-frog time scheme. This kind of method provides a dissipation-free approximation of propagation equations and allows for the accurate estimation of aeroacoustic energy variation, which is not possible with numerical methods (finite volumes, discontinuous Galerkin, spectral elements) based on upwind numerical fluxes. The fact that centered numerical fluxes in discontinuous Galerkin time-domain methods can lead to discretization methods inducing no numerical dissipation is quite well known. This was for example first 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 uniform supporting flow, at the continuous level (i.e. before space discretization), some quadratic energy verifies a balance equation without source term. This means energy is conserved (up to boundaries); 2. in this “uniform supporting flow” case, we are able to prove that our DG method (with leap-frog time-scheme and centered fluxes) 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 claim that we “control energy variations in the uniform case”; 3. accordingly, for a non-uniform supporting flow, at the continuous level (i.e. before space discretization), we use the well-known symmetrization of nonlinear Euler equations 17 to derive an aeroacoustic energy which verifies some balance equation with source term. Because of this unsigned source term, aeroacoustic waves can be damped or excited by the supporting flow. 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 uniform supporting flow” 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 term is also verified. 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 them artificially (like methods based on upwind fluxes, which damp instabilities, in an uncontrolled way though); 5. we show finally that, by introducing an artificial 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 uniform or non-uniform supporting flow, 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 different configurations, obtained with MPI-based Fortran parallel implementations 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 flow. Therefore, we linearize the three-dimensional Euler equations around a given steady flow and only take into account first-order perturbation terms. For a perfect inviscid gas, Euler equations read:

∂t

ρ ρu ρv ρw e

+ ∂x

ρu ρu2 + p ρuv ρuw (e + p)u

+ ∂y

ρv ρuv ρv 2 + p ρvw (e + p)v

+ ∂z

ρw ρuw ρvw ρw2 + p (e + p)w

= 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 − 12 ρk~v k2 ), where γ is a fixed constant (γ > 1). 2.1. Linearization around a uniform flow A first step towards aeroacoustic wave propagation consists in linearizing Euler equations around any uniform flow defined by some constant physical fields (ρ 0 , ~v0 , p0 ), therefore being a steady-state solution of Euler equations. Infinitely small perturbations (δρ, δ~v , δp) of this flow verify the following linear hyperbolic system of conservation laws: ~ + A x ∂x W ~ + A y ∂y W ~ + A z ∂z W ~ = ~0, ∂t W

(2.2)

t

~ = (δp − c2 δρ, t ρ0 c0 δ~v , δp), c0 corresponds to the uniform sound speed given by where W 0 2 ~ has been chosen such that the constant matrices A x , Ay , and c0 = γp0 /ρ0 . The variable W Az are symmetric. They are given by: u0 0 0 0 0 v0 0 0 0 0 w0 0 0 0 0 0 u0 0 0 c 0 0 v0 0 0 0 0 w0 0 0 0 Ax = 0 0 u0 0 0 , Ay=0 0 v0 0 c0 , Az= 0 0 w0 0 0 0 0 0 u0 0 0 0 0 v0 0 0 0 0 w 0 c0 0 c0 0 0 u 0 0 0 c 0 0 v0 0 0 0 c 0 w0 In the particular case of the linearization around a uniform flow, one can define a volumic ~ 2 (a very simple expression in function of W), ~ which verifies aeroacoustic energy E = 21 kWk ~ = 0, where the energy flux F ~ is given by the balance equation ∂t E + divF

1t ~ ~ s ∈ {x, y, z}. WAs W, 2 Thus, away from boundaries, the aeroacoustic energy E is exactly conserved. Fs =

2.2. Linearization around a non-uniform flow A more interesting framework consists in linearizing Euler equations around a nonuniform steady-state solution. In that case, the steady flow is defined by smoothly varying physical quantities (ρ0 , ~v0 , p0 ). Linearizing straightforwardly Euler equations (2.1) yields: ~ + ∂x A0x W ~ + ∂y A0y W ~ + ∂z A0z W ~ = 0, ∂t W (2.3)

~ now denotes the perturbations of conservative variables (i.e. W ~ T = (δρ, ρ0 δ~v + where W ~v0 δρ, δp/(γ − 1) + ρ0~v0 .δ~v + k~v0 k2 δρ/2)) and the space-varying matrices A 0x , A0y , and A0z are given in function of γ˜ = γ − 1, α0 = c20 /˜ γ + k~v0 k2 /2, β0 = (γ − 2)kV0 k2 /2 − c20 /˜ γ , and 3 the canonical basis (~ex , ~ey , ~ez ) of R by t 0 ~es 0 A0s = γ2˜ k~v0 k2~es − (~v0 .~es )~v0 (~v0 .~es )I3 − γ˜~es t~v0 + ~v0 t~es γ˜~es , s ∈ {x, y, z}. (2.4) t t β0 (~v0 .~es ) α0 ~es − γ˜ (~v0 .~es ) ~v0 γ(~v0 .~es )

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

In this equation, the matrices A0x , A0y , and A0z are not symmetric anymore, and it is very difficult to deduce any aeroacoustic energy balance equation. Therefore, we consider other acoustic variables, derived from the quite classical symmetrization of Euler equations. Assuming the flow is smooth enough, the change of variables (ρ, ρ~v , e) → γ γ e˜ γ v , − ρ˜ (− p + γ + 1 − ln ρpγ , ρ˜ p~ p ) transforms Euler equations (2.1) into a symmetric system of conservation laws (i.e. Jacobians of fluxes 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 ~ + ∂x A ˜ 0V ~ + ∂y A ˜ 0V ~ + ∂z A ˜ 0V ~ = 0, A00 ∂t V (2.5) x y z ~ is given in function of variables W ~ by V ~ = A0 −1 W ~ and where V 0 t 1 ~v0 α0 − c20 /γ ρ0 ˜0 c20 A00 = ~v0 v0 t~v0 α0~v0 ; As = A0s A00 , s ∈ {x, y, z}. (2.6) γ I3 + ~ γ˜ α0 − c20 /γ α0 t~v0 α20 − c40 / (γ˜ γ)

A00 clearly is symmetric and it can be proved that it is definite positive (and then not ~ by A0 V ~ in Eq. 2.3 (and by singular). Eq. 2.5 can also be obtained simply by replacing W 0 0 ˜s noting that ∂t A0 = 0). Finally, the reader can also check that the symmetric matrices A (s ∈ {x, y, z}) are given by t 0 ~es (~v0 .~es ) ˜ 0 = (~v0 .~es ) A0 + p0 ~es (~v0 .~es )~v0 + α0~es . ~es t~v0 + ~v0 t~es A s 0 γ˜ t t (~v0 .~es ) (~v0 .~es ) ~v0 + α0 ~es 2α0 (~v0 .~es )

t ~ 0 −1 W ~ ≡ 1 t VA ~ 0V ~ verifies Then, the volumic aeroacoustic energy E defined by E = 21 WA 0 0 2 the following balance equation with source term: ( t ~A ˜ 0 V, ~ s ∈ {x, y, z}. Fs = V sh i ~ ∂t E + divF = S, with (2.7) t ~ ~ ∂x (A ˜ 0 ) + ∂ y (A ˜ 0 ) + ∂ z (A ˜ 0 ) V. S = − 21 V x y z

Thus the aeroacoustic energy is not conserved and the variations in the steady flow considered can damp or amplify aeroacoustic waves, unless the source term vanishes (which is the case for a uniform flow for example). In the sequel, we shall mainly discretize the conservative 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 flow 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 field artificial boundary ∂Ω∞ ). Ω is partitioned into a finite number of polyhedra (each one having a finite number of faces). For each polyhedron T i , called ”control volume” or ”cell”, Vi 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 T T F int = F/∂Ω, physical faces F phys = F ∂Ωphys and absorbing faces F abs = F ∂Ω∞ . For T each internal face aik = Ti Tk , we denote by Sik the measure of aik and by ~nik the unitary normal, oriented from Ti towards Tk . The same definitions are extended to boundary faces, the index k corresponding to a fictitious 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 P face in common). We also define the perimeter P i of Ti by Pi = k∈Vi Sik . We recall the P nik = 0. following geometrical property for all control volumes: k∈Vi Sik ~ Following the general principle of discontinuous Galerkin finite element methods, the unknown field inside each control volume is seeked for as a linear combination of local basis vector fields ϕ ~ ij , 1 ≤ j ≤ di (generating the local space Pi ) and the approximate field is allowed to be fully discontinuous across element boundaries. Thus, a numerical flux function has to be defined to approximate fluxes at control volumes interfaces, where the approximate solution is discontinuous. This context is quite general. Actual implementations of the method have been considered only on tetrahedral meshes, where control volumes are the tetrahedra themselves. We shall only consider constant (P0 ) or linear (P1 ) approximations inside tetrahedra. 3.1. Time and space discretizations We only consider here the most general case of aeroacoustic wave propagation in a nonuniform steady flow. Also, in order to limit the amount of computations, we restrict our study to piecewise constant matrices A 0s (s ∈ {x, y, z}) given in Eq. (2.4). For each control volume Ti , for s ∈ {x, y, z}, we denote by Ais an approximate for the average value of A 0s over Ti . Dot-multiplying Eq. (2.3) by any given vector field ϕ ~ , integrating over T i and integrating by parts yields Z

Z Z X X ~ ∂W t ~ − ~ . ϕ ~· ∂s ϕ ~ A0s W = ϕ ~· ns A0s W ∂t Ti Ti ∂Ti s∈{x,y,z}

(3.8)

s∈{x,y,z}

~ by the approximate field W ~ i and Inside volume integrals over Ti , we replace the field W ~ the matrices A0s by their respective average values A is . For boundary integrals over ∂Ti , W

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

is discontinuous, and we define totally centered numerical fluxes, i.e.: h i ( ~ i + Pk W ~ k , ~ ∀i, ∀k ∈ Vi , (nikx A0x + niky A0y + nikz A0z )W ' 21 Piik W ik |aik

with Piik = nikx Aix + niky Aiy + nikz Aiz ,

Pkik = nikx Akx + niky Aky + nikz Akz .

(3.9)

Concerning the time discretization, we use a three-level leap-frog scheme. The unknowns ~ i are approximated at integer time-stations t n = n∆t. Assuming we dispose of W ~ n−1 W i n+1 ~ n , the unknowns W ~ and W are seeked for in Pi such that, ∀~ ϕ ∈ Pi , i i Z

Ti

ϕ ~·

~ n+1 − W ~ n−1 Z W i i = 2∆t Ti

X

t

~ in − ∂s ϕ ~ Ais W

s∈{x,y,z}

XZ

ϕ ~·

aik

k∈Vi

~ n + Pk W ~ n Piik W i ik k . 2

(3.10)

Again, the time scheme above is almost explicit. Each time step only requires the inversion of local symmetric positive definite matrices of size (d i × di ). In the particular case where Pi is a complete linear (P1 ) 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 flow and the acoustic perturbations. This means that we assume that for any slip boundary face a ik belonging to the control volume Ti , the steady solution of Euler equations verifies a slip condition at the discrete level, i.e. ~nik · ~vi0 = 0. For the acoustic perturbations, we use a mirror fic~ k in the computation of the boundary flux given in Eq. (3.10). We take titious state W δρk = δρi , δpk = δpi , and δ~vk = δ~vi − 2(~nik · ~vi )~nik (which implies (δ~vk − δ~vi ) × ~nik = 0 and δ~vk .~nik = −δ~vi .~nik ). For an absorbing boundary face aik , upwinding is used to select outgoing waves only. Before discretization in time, classical upwinding leads to a boundary flux F ik given by + i ~ Fik = Pik Wi , where for any diagonalizable matrix Q = S −1 DS with D diagonal, Q+ = (Q + |Q|)/2 and terms of the diagonal matrix |D| are the moduli of the eigenvalues. This ~ i . However, for this intuitive numerical flux, it is ~ k = |Pi |W general idea leads to Pkik W ik very difficult to prove that the resulting time-scheme is stable and that energy is actually sent in the exterior domain. We then consider the numerical flux based on the following √ √ i √ i −1 i √ i √ i −1 W ~ n−1 +W ~ n+1 n k i i ~ , where Ai is the positive fictitious state: P W = A A P A A ik

k

0

0

ik

0

0

2

0

square root of the symmetric definite positive matrix A i0 . Indeed, this expression derives from the intuitive upwind flux 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 √ i √ i −1 i √ i √ i −1 n−1 n 22 k ~ = A A P A A W ~ available . It takes the form Pik W . 0 0 0 0 i k ik

Bernacki, Piperno

3.3. Energy balance and stability In order to investigate stability, we define a discrete aeroacoustic energy F n by: Z t 1X n ~ n Ai −1 W ~ n + tW ~ n+1 Ai −1 W ~ n−1 W F = i 0 i 0 i i 2 T i i Z t ∆t X −1 ~ n−1 i −1 ~ n−1 n+1 ~ + , M A W + W Ai0 W ik 0 i i i 4 abs aik aik ∈F

√ √ −1 i √ i −1 √ i √ √ −1 √ √ ˜ A A . One can show that the with Mik = Ai0 Ai0 Piik Ai0 Ai0 ≡ Ai0 Ai0 P 0 0 ik matrices Mik 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 − Fn ): Z t n k ∆t X n+1 n ~ (P ˜ −P ˜ i )V ~ n+1 + t V ~ n+1 (P ˜k − P ˜ i )V ~n V F −F = − i ik ik ik ik k i k 2 aik aik ∈F int Z t ∆t X ~ n−1 + V ~ n+1 Mik V ~ n−1 + V ~ n+1 . − V (3.11) i i i i 4 abs aik aik ∈F

~ n ≡ Ai −1 W ~ n , ∀i, ∀n. The first term is a where we have used the auxiliary variables V 0 i i 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 influence on the global energy balance. In order to prove stability, we show in Annex 2 that F n is a quadratic positive definite ~ n ) under some CFL-like sufficient stability condition ~ n−1 , W form of numerical unknowns (W i i on the time-step ∆t: 2Vi ∀i, ∀k ∈ Vi , ∆t (2λi αi + βik ρik ) < , (3.12) Pi where αi and βik are dimensionless regularity coefficients depending of basis functions and element aspectratio, λi = |ui0 | + |v0i| + |w0i | +3ci0 , and ρik = |~v0i ·~nik | + ci0 fora boundary face 2 2 −1 −1 for an internal , |~v0k · ~nik | + ck0 ρ Ai0 Ak0 and ρ2ik = sup |~v0i · ~nik | + ci0 ) ρ Ak0 Ai0 face (ρ here denotes the spectral radius of a matrix). In the case of a uniform flow, we have P iik = Pkik . Thus the aeroacoustic energy is non-increasing (and exactly conserved if no absorbing boundary is present, which shows the scheme is genuinely non-diffusive) and the scheme is stable under a CFL-type stability condition depending on the size of elements and sup i (k~vi k + c0i ). 4. Numerical results We dispose of a three-dimensional parallel implementation of the DGTD method presented in the previous section. Any subsonic steady flow can be considered, even with strong flow gradients. However, the flow, given as the output of a non-linear Euler equations solver, has to be post-processed: average of the flow 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 testcases in two and three space dimensions, in order to validate the method on benchmark problems, test the method on complex flows and configurations, and finally evaluate the performance of the parallel Fortran 77 implementation, based on the MPICH implementation 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 efficiency. 4.1. Acoustic perturbations in a three-dimensional uniform flow This quite classical test-case23,24 is aimed at validating the DGTD method proposed by mixing acoustic, entropic and vorticity perturbations, in a uniform horizontal flow (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 M1 M2

# vertices 68,921 531,441

# tetrahedra 384,000 3,072,000

# absorbing faces 19,200 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 finer mesh M2. Artificial reflections 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 flow is uniform, the total acoustic energy is first exactly conserved, and then decreasing, once perturbations have reached the absorbing boundaries (at t = 15). Computation times are given in Table 2. Parallel efficiency reaches a satisfying level and total CPU time are really acceptable. 4.2. Acoustic perturbations in a three-dimensional subsonic Couette flow This second test-case is more complex and is classically used to evaluate the merits of a numerical method for non-uniform shear flows 25 . The steady flow considered is horizontal, with ~v0 = 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 flow. 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.

3.5 3

F

2.5 2

1.5 1 0.5

0

5

10 15 20 25 30 35 40 45 50 55 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 efficiency. Mesh M1 M2 -

Np 1 4 8 16 8 16

CPU 33381 s 8167 s 4286 s 2195 s 34012 s 17203 s

REAL 33447 s 8412 s 4465 s 2381 s 35651 s 18760 s

% CPU 100% 97% 96% 92% 95% 92%

S(Np ) 1 4 7.5 15.2 1 1.98

waveguide, the aim of this test-case is to validate the accurate propagation and convection of the waves by the shear flow, 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 flow test-case. Characteristics of meshes. Mesh M1 M2

# vertices 73,629 522,801

# tetrahedra 405,600 3,000,000

# absorbing faces 23,504 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 flow : in the numerical results, artificial reflections from absorbing boundaries have a relative amplitude of 1.5% for the coarse mesh M1 (respectively 1% for the finer 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 flow (thus the bottom part of this perturbation travels slowly with the supporting flow, whose velocity vanishes in the plane z = 0). The time evolution of the global acoustic energy is given on Fig. 4. As expected, the energy is not conserved anymore because the non-uniform flow provides energy to the acoustic perturbations. The energy decreases when perturbations reach the absorbing boundaries. A part of the energy corresponding to the bottom part of the entropic perturbation remains in the domain. Computation times are given in Table 4 and reveal similar parallel efficiency 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 flows. 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 flow 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.08 0.07 0.06

Energy

0.05 0.04 0.03 0.02 0.01 0

0

1

2

3 Time (s)

4

5

Fig. 4. Couette flow test-case. Time evolution of the acoustic energy.

Table 4. Couette flow test-case. Parallel CPU times and efficiency. Mesh M1 M2 -

Np 4 8 16 8 16

CPU 12830 s 6428 s 3230 s 46359 s 23412

REAL 13008 s 6590 s 3373 s 47743 s 24959 s

% CPU 99% 98% 96% 97% 94%

S(Np ) 1 2 3.9 1 1.98

Bernacki, Piperno

acoustic energy in the flow. 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 flow 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, p0 = 1/γ. In this flow, u0 goes through an inflection point at y = 0.25 (and k~v 0 k = 0.5 on the inflection point). In this particular flow, the acoustic energy is amplified along the inflection line and at the same time, perturbations are convected along this line, which leads to an instability. The behavior is completely different from the one obtained with a shear flow of the form ~v0 = (y + 0.25)~ex , 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 flow with an inflection point. Different means are under current investigation to control this kind of instability.

Fig. 5. Contours of kδ~v k at t=0.257 for the linear (left) and inflection (right) shear flows with source term

4.4. Two-dimensional instabilities past a NACA airfoil profile This test-case is aimed at evaluating the occurrence of Kevin-Helmholtz-type instabilities for steady-state flows with sharp gradients. We have considered a test-case past a NACAtype airfoil profile (courtesy of ONERA 19 ). The steady-state solution at M∞ = 0.5 obtained by our finite-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 profile, and our absorbing boundary conditions in the far field boundary. For this testcase, 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δ~v k are shown on Fig. 7. Some numerical treatments are possible to overcome these difficulties and get rid of instabilities. One of them4 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δ~v k are shown on the right part of Fig. 8 and can be compared to classical acoustic scattering (uniform still flow). The comparison of contours near the trailing edge for acoustics and aeroacoustics is shown on Fig. 9, where the shear flow at the trailing edge has a visible influence 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 flow past the airfoil profile at M∞ = 0.5

Fig. 7. Contours of kδ~v k at t = 1.57ms for the original aeroacoustic model.

Bernacki, Piperno

propagation of acoustics (and generation of vortices).

Fig. 8. Contours of kδ~v k for the stabilized aeroacoustic model (left) and classical acoustics (right).

Fig. 9. Tailing edge zoom – kδ~v k 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 flow past a sphere and the propagation of waves emitted by an acoustic source (Gaussian profile, periodic source term with period T = 0.2s. The subsonic flow with M ∞ = 0.5 is obtained by a three-dimensional MUSCL-extended finite-volume solver of three-dimensional Euler equations on tetrahedral grids. The steady-state solution presents Mach numbers between 0.002 and 0.8 and contours 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 field 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 324,471

# elements 1,870,288

# abs. faces 10,092

# slip faces 46,080

Np 16

CPU 50h4mn

REAL 51h18mn

% CPU 97,6%

resulting propagated waves for this aeroacoustic case in a steady flow past a sphere (after 10 and 25 periods). The numerical results (scattering and reflection 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 configurations. 5. Conclusion In this paper, we have presented a time domain Discontinuous Galerkin dissipationfree 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 flow, we have proved, using the well-known symmetrization of Euler equations, that the aeroacoustic energy satisfies a balance equation with source term at the continuous level, and that our numerical method satisfies an equivalent balance equation at the discrete level. This shows that the method is genuinely dissipation-free, which is confirmed by numerical results, including test-cases where Kelvin-Helmholtz instabilities appear. Further works can concern many different 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 flows 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 flow (currently, it is only P0 ). 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 i0 is symmetric, we have

F

n+1

−F

n

=

1X 2 i

∆t + 4

Z

t

Ti

~ n+1 Ai −1 W ~ n+1 − W ~ n−1 ~ n Ai −1 W ~ n+2 − W ~ n + tW W 0 i 0 i i i i i

X

aik ∈F abs

Z

aik

t

~ n+2 + W ~ n ~ n Mik Ai −1 W W 0 i i i

Ai0

−1

t

−1 ~ n−1 Ai0 W i

−

−1 Mik Ai0

~ n+1 + W ~ n−1 W i i

.

The variational form of the scheme (3.10) allows to consider any field in P i , thus for example

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

~ n−1 or W ~ n themselves. We obtain that Fn+1 = Fn + ∆t(T1 − 1 T2 − 1 T3 ) with the fields W i i 2 2 XZ

T1 =

i

Ti s∈{x,y,z}

X Z

T2 =

aik ∈F int

aik ∈F int

T3 =

X

aik ∈F slip

+

" X Z t

aik ∈F abs

aik

Z

t

aik

t

aik

X Z

+

t

X

t ˜iV ~n ~ n+1 A ˜iV ~ n+1 + ∂s V ~n A ∂s V s i s i i i

~ in P ˜i V ~ n+1 + P ˜k V ~ n+1 + t V ~ n+1 + P ~ n+1 ˜k V ˜i V ~n P V ik i ik k ki k ki i k t

aik

t ~ n+1 P ˜k V ~ n ˜i ~ n ~ n+1 P ˜i V ~ n ˜k ~ n V ki k + Pki Vi ik i + Pik Vk + Vk i

˜ i Hik V ~n ˜i V ~n+P ~ n+1 P ˜i V ~ n+1 + P ˜ i Hik V ~ n+1 + t V ~n P V ik i ik i ik i ik i i i

! !# ~ n+ V ~ n+2 ~ n−1 + V ~ n+1 V V t i i i ˜i V ˜i V ~ n+1 + Mik i ~ n+1 P ~n P + V ik i ik i + Mik i 2 2 Z t n 1 X ~ n−1 Mik V ~ n−1 + V ~ n+1 . ~ Mik V ~n+V ~ n+2 − t V − V i i i i i i 2 abs aik

~ in V

aik ∈F

˜ is symmetric, and since P ˜ i = −P ˜ i and P ˜ k = −P ˜ k , we obtain that Since P ik ki ik ki X Z

T2 = 2

aik

aik ∈F int

+

t

X Z

aik

aik ∈F int

T2 + T 3 = 2

XXZ i

+

k∈Vi

aik ∈F int

+

1 2

t aik

X Z X

t

aik ∈F abs

Z

~ n+1 , and ˜i V ˜k − P ~ n+1 + t V ~n P ˜i V ˜k − P ~n P V i ik ik k ik ik i k

~ nP ˜ i ~ n+1 + V i ik Vi

X

aik ∈F slip t

aik

~ nP ˜ i ~ n+1 + t V ~ nP ˜ k ~ n+1 V i ik Vi k ki Vk

Z

t

aik

~ nP ˜i ~ n+1 + t V ~ n+1 P ˜ i Hik V ~n V i ik Hik Vi ik i i

~ n+1 ~ n+1 + t V ˜i V ˜i V ~n P ˜k − P ~n P ˜k − P V i ik ik ik ik k i k t

aik

~ n−1 + V ~ n+1 Mik V ~ n−1 + V ~ n+1 . V i i i i

In the expressions above, the matrix H ik denotes the matrix used for slip boundary con~ n is given by W ~ n = Hik W ~ n (Hik does not ditions, in the sense that the fictitious state W i k k change the density and energy perturbations, and operates an orthogonal symmetry on the velocity perturbation). One can also show that H 2ik = I and that Hik commutes with Ai0 .

Bernacki, Piperno

We finally get that Fn+1 = Fn + ∆t(T4 + T5 + T6 ) with XZ X Z X t t t n i ~ P ˜ ~ n+1 ˜iV ~n − ˜iV ~ n+1 + ∂s V ~ n+1 A ~n A V T4 = ∂s V i ik Vi s i s i i i Ti s∈{x,y,z}

i

T5 = −

1 2

T6 = −

1 4

−

1 2

X Z

t

aik

aik ∈F ref

~ in P ˜ i Hik + V ik

Z

t

X Z

t

X

aik ∈F abs

aik ∈F int

aik

aik

k∈Vi

t

˜ i Hik P ik

aik

~ n+1 V i

~ n−1 + V ~ n+1 Mik V ~ n−1 + V ~ n+1 V i i i i

˜k − P ˜i V ˜k − P ˜i V ~ n+1 ~ n+1 + t V ~n P ~ in P V ik ik k ik ik i k

˜ s are symmetric and constant over control volumes, we get by integratSince the matrices A ˜ i Hik is skew-symmetric, we get T5 = 0. ing by parts that T4 = 0. Also, since the matrix P ik ~ variables, we get the result announced. Using the expression of Mik and using the W

Annex 2 ~ n−1 and W ~ n We shall prove that Fn is a positive definite quadratic form of unknowns W i i under the condition given in Eq. (3.12). Let us introduce the following definitions. R ~ ∈ Pi , we denote by kXk ~ T = ( kXk ~ 2 )1/2 the L2 -norm of X ~ over Ti . Definition 1 ∀X i Ti R 2 1/2 ~ ~ a =( . Finally, the matrix Ai0 For any interface aik ∈ F, we also define kXk ik aik kXk ) √ ~ T. ~ T = k Ai −1 Xk being positive definite, we define a second norm by k| Xk| i i 0 Definition 2 We assume some regularity for the basis functions ϕ ~ ij , 1 ≤ j ≤ di . More precisely, we assume that, for each control volume T i , there exist dimensionless positive constants αi and βik (for each k ∈ Vi ) such that ( αi Pi ~ ~ ~ ∈ Pi , ∀s ∈ {x, y, z}, k∂s XkTi ≤ Vi kXkTi , ∀X (5.13) ~ 2. ~ 2 ≤ βik k~nik k kXk ∀aik ∈ F, kXk aik

Vi

Ti

√ ~ n = Ai −1 W ~ n . We introduce Definition 3 We choose here some notations. We define Z 0 i i P i i i i i λ = s ρ As (and we recall ρ As = |~es .~v0 | + c0 ). For boundary faces aik , we recall that some neighboring fictitious control volume T k is imagined, and we take by convention: ~ n k|T , βki = βik , Vk = Vi . Finally, we define ρik for ~ n k|T = k|W ~ n kT , k|W ~ n kT = k W kW i i i k i k k k any face aik by: 2 2 −1 −1 , |~v0k · ~nik | + ck0 ρ Ai0 Ak0 aik ∈ F int , ρ2ik = sup |~v0i · ~nik | + ci0 ρ Ak0 Ai0 aik ∈ F slip , ρik = ρ Piik = |~v0i · ~nik | + ci0 aik ∈ F abs , ρik = 0,

A dissipation-free DGTD method for the three-dimensional linearized Euler equations

~ n−1 and W ~ n In order to prove that Fn is a positive definite quadratic form of unknowns W i i under the condition (3.12), we first show the following result: i 1 ~ n 2 ~ n−1 k|2 − αi Pi λ ∆t k|W ~ n k|T k|W ~ n−1 k|T k|Wi k|Ti + k|W Fni ≥ Ti i i i i i 2 Vi ∆t X ρik βik k~nik k ~ n−1 2 ρki βki k~nki k ~ n 2 − (5.14) k|Wi k|Ti + k|Wk k|Tk . 4 Vi Vk k∈Vi

We can show easily that the above lemma gives the result. Using the above definitions and P ρki ∆t ~ n 2 assumptions, we have Fni ≥ k∈Vi k~nik k(Xik − βki4V k|Wk k|Tk ) with k i 1 βik ρik ∆t 1 ~ n 2 ~ n k|T k|W ~ n−1 k|T ~ n−1 k|2 − αi λ ∆t k|W k|Wi k|Ti + − k|W Xik = i Ti i i i i 2Pi 2Pi 4Vi Vi 1 αi λi ∆t 1 βik ρik ∆t αi λi ∆t ~ n k|2 + ~ n−1 k|2 − − − ≥ k|W k|W i Ti Ti i 2Pi 2Vi 2Pi 4Vi 2Vi X X k~nik kZik with k~nik kYik + Finally, Fn ≥ aik ∈F int

Yik

Zik

aik ∈F slip

S

F abs

1 βik ρik ∆t αi λi ∆t ~ n 2 ~ n−1 k|2 − − k|Wi k|Ti + k|W = Ti i 2Pi 4Vi 2Vi βki ρki ∆t αk λk ∆t ~ n 2 1 ~ n−1 k|2 , − − k|Wk k|Tk + k|W + Tk k 2Pk 4Vk 2Vk βik ρik ∆t αi λi ∆t ~ n 2 1 ~ n−1 k|2 . = − − k|Wi k|Ti + k|W Ti i 2Pi 4Vi 2Vi

Therefore, the energy Fn is a positive definite quadratic form of unknowns under the sufficient CFL-type stability condition proposed in Eq. (3.12). End ofthe proof. We now give the proof for the preliminary result (5.14). We first have n−1 2 1 n 2 n n n ~ ~ Fi = 2 k|Wi k|Ti + k|Wi k|Ti + ∆t 2 (Xi + Yi ) with Z X t 1 −1 ~ n−1 i −1 ~ n−1 n+1 n ~ M A W + W Ai0 W Xi = ik 0 i i i 2 aik aik ∈F abs ∩∂Ti Z X t n−1 1 ~ ~ n−1 + V ~ n+1 = V M V ik i i i 2 aik aik ∈F abs ∩∂Ti Z t 1 n ~ n−1 Ai −1 W ~ n+1 − W ~ n−1 W Yi = 0 i i i ∆t Ti Z X t XZ t n−1 i n ˜ V ~ − ~ ˜i V ~n+P ˜k V ~n ~ n−1 P = 2 A ∂s V V i

Ti s∈{x,y,z}

=

Z

X

Ti s∈{x,y,z}

s

i

k∈Vi

t

aik

i

ik

i

ik

XZ t n−1 i n n ˜ ∂s V ˜iV ~ ~ ~ ~ n−1 A − A − V ∂s V s i s i i i k∈Vi

aik

t

k

˜k V ~n ~ n−1 P V ik k i

Bernacki, Piperno

Using Hik and eliminating terms corresponding to absorbing boundaries leads to Z X t t n−1 i i n n n1 n n ˜ V ˜ ∂s V ~ ~ ~i ~i A A Xi + Y i = ∂s V s i − Vi s Ti s∈{x,y,z}

−

X Z

aik ∈F int

t aik

~ n−1 P ˜k V ~n V ik k − i

X

aik ∈F slip

Z

t aik

~ n. ~ n−1 P ˜ i Hik V V i ik i

Then we deduce that Xni + Yni = T1 + T2 + T3 with Z √ X t √ √ √ ~ n−1 Ai −1 A ˜ i Ai −1 Z ~ n − tZ ~ n−1 Ai −1 A ˜ i Ai −1 ∂s Z ~n , ∂s Z T1 = 0 s 0 i 0 s 0 i i i Ti s∈{x,y,z}

T2 = − T3 = −

Z

X

aik ∈F int ∩∂Ti

X

aik ∈F slip ∩∂Ti

t aik

Z

aik

√ √ √ √ ˜ k Ak −1 Z ~ n, ~ n−1 Ai −1 Ak Ak −1 P Z 0 0 0 ik 0 k i

t

√ √ √ √ ˜ i Ai −1 Ai Hik Ai −1 Z ~ ni . ~ n−1 Ai0 −1 P Z 0 0 0 ik i

√ −1 √ −1 √ −1 i √ i −1 ˜ A ˜ i Ai = ρ Ais , then we have The matrix Ai0 A is symmetric and ρ Ai0 A s s 0 0

2λi αi Pi ~ n−1 ~ n kT = 2λi αi Pi k|W ~ n−1 k|T k|W ~ n k|T . k Z i k Ti k Z i i i i i i Vi Vi √ −1 √ −1 √ −1 k √ k −1 ˜ A ˜ k Ak = ρ Pkik , then is also symmetric and ρ Ak0 P The matrix Ak0 P 0 0 ik ik |T1 | ≤

|T2 | ≤ ≤

ρ

Pkik

aik ∈F int ∩∂Ti

r √ −1 √ k ~ n ka ~ n−1 ka kZ ρ Ak0 Ai0 A0 kZ i ik ik i

X

ρ

Pkik

r −1 ~ in k|a . ~ n−1 k|a k|W k|W ρ Ak0 Ai0 ik ik i

X

aik

∈F int ∩∂T

i

Since the matrices Ai0 and Hik commute and H2ik = Id, we have X X n−1 ~ n k|a . ~ n−1 k|a k|W ~ n ka = ~ ρ Piik k|W kaik kZ |T3 | ≤ ρ Piik kZ i i ik ik ik i i aik ∈F ref ∩∂Ti

aik ∈F ref ∩∂Ti

1 2 + b2 , rewriting Fn ≥ Using the above definitions and inequalities of the form ab ≤ a i 2 n−1 2 ∆t 1 2 n ~ ~ k|Ti − 2 (|T1 | + T2 | + T3 |) leads to the intermediate result (5.14). 2 k|Wi k|Ti + k|Wi 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. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25.

M. Lesieur and O. Metais, Annu. Rev. Fluid Mech., 28, pp.45-82 (1996). C.K.W. Tam and J.C. Webb, J. Comput. Phys., 107, pp.262–281 (1993). C. Bogey, C. Bailly and D. Juv´e, AIAA J., 40(2), pp.235-243 (2002). C.K.W. Tam, AIAA J., 33(10), pp.1788-1796 (1995). M.J. Lighthill, Proc.Roy. Soc. London, 211, Ser.A, 1107, pp.564-587 (1952) 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). C. Bogey, Calcul direct du bruit a´erodynamique et validation de mod`eles acoustiques hybrides, Ph.D. thesis, Ecole centrale de Lyon, 2000. L. Fezoui, S. Lanteri, S. Lohrengel and S. Piperno, RAIRO M2AN, 39, (6)pp.1149-1176 (2005) F.Q. Hu, J. Comput. Phys., 129, pp.201-219 (1996) F.Q. Hu, J. Comput. Phys., 173, pp.455-480 (2001) F.Q. Hu, M.Y. Hussaini and P. Rasetarinera, J. Comput. Phys., 151, pp.921-946 (1999) D. Kr¨ oner, Math. Comput., 195(57) (1991) E. Becache, S.D. Fauqueux and P. Joly INRIA Research Report RR-4304 (INRIA, 2001). C.K.W. Tam and Z. Dong, J. Comput. Phys., 4(2), pp.175-201 (1996). S. Piperno and L. Fezoui, INRIA Research Report RR-4733 (INRIA, 2003) A. Harten, ICASE Research report, (81-34) (1981). J. Hesthaven and T. Warburton, J. Comput. Phys. 181 (1) 186–221 (2002). P. Delorme and C. Peyret, Galerkin discontinous method for computational aeroacoustics, Private communication. M. Bernacki and S. Piperno, INRIA Research Report RR-4932 (INRIA, 2003). M. Bernacki and S. Piperno, INRIA Research Report RR-4699 (INRIA, 2003). M. Bernacki, S. Lanteri, and S. Piperno, to appear in J. Comput. Acoust.. C.K.W Tam and J.C. Hardin, Second computational aeroacoustics workshop on benchmark problems, NASA Ohio Aerospace Institute (1997). C.K.W Tam, J.C. Hardin and J.R. Ristorcelli, Workshop on benchmark problems in computational aeroacoustics, NASA Ohio aerospace institute (1995). Fourth computational aeroacoustics workshop on benchmark problems NASA, Ohio Aerospace Institute (2003).