A. Crivellini, F. Bassi 1

SVILUPPO DI UN METODO AGLI ELEMENTI FINITIDISCONTINUI DI GALERKIN PER LO STUDIO DELLA

PROPAGAZIONE ACUSTICA

Andrea Crivellini, Francesco Bassi

Dipartimento di Energetica, Universita degli Studi di Ancona

SOMMARIO

In questo lavoro viene utilizzato un metodo agli elementi finiti discontinui di Galerkin per lostudio della propagazione acustica descritta dalle equazioni di Eulero linearizzate attorno aduno stato medio non uniforme.Il metodo, dotato di elevato ordine di accuratezza e adatto a griglie non strutturate, ha ilgrande vantaggio di essere insensibile alle caratteristiche della griglia. Permette quindi l’utilizzodi reticoli di calcolo con forti addensamenti e con brusche variazioni di dimensioni degli elementisenza che l’accuratezza della soluzione ne risenta.La soluzione avanza nel tempo con un metodo di Runge–Kutta esplicito a cinque passi, accuratoal quarto ordine, che necessita solo di due registri per la memorizzazione.Una tecnica di buffer layer e utilizzata per troncare il dominio in un dominio di calcolo finitosenza introdurre riflessioni spurie dovute alla presenza dei contorni.Il codice di calcolo e stato applicato alla soluzione di un caso test bidimensionale della Compu-tational Aeroacoustics (CAA) caratterizzato da un flusso medio non uniforme, rappresentanteun getto supersonico, e dall’introduzione nelle equazioni di un termine sorgente che eccita lapropagazione ondosa.I risultati dimostrano la capacita del metodo di calcolare la corretta propagazione acustica anchedelle perturbazioni il cui valore risulta essere di molti ordini di grandezza inferiore a quello dellegrandezze medie.

1 EQUAZIONI DI EULERO LINEARIZZATE

Le equazioni di Navier-Stokes, in termini delle variabili primitive q = (ρ, p, u, v)T e per un gasideale, possono essere scritte come:

∂q

∂t+ A

∂q

∂x+ B

∂q

∂y= d . (1)

Il termine d contiene i termini dovuti alla viscosita, mentre A e B sono le due seguenti matrici:

A =

u 0 ρ 00 u γp 00 1/ρ u 00 0 0 u

B =

v 0 0 ρ0 v 0 γp0 0 v 00 1/ρ 0 v

. (2)

Scomponendo il vettore q nella somma tra la soluzione q di riferimento e la deviazione q′ daquesta,

2 56 Congresso Nazionale ATI

q (x, y, t) = q (x, y) + q′ (x, y, t) , (3)

trascurando i termini in q′ di ordine superiore al primo, si ottiene la seguente linearizzazioneattorno allo stato q:

∂q′

∂t+ A

∂q′

∂x+ B

∂q′

∂y+ Cq′ − d′ = −∂q

∂t− A

∂q

∂x− B

∂q

∂y+ d . (4)

Le matrici A e B sono le (2) calcolate in base allo stato di riferimento e la matrice C e uguale a

C =

ux + vy 0 ρx ρy

0 γ (ux + uy) px py

(uux/ρ + vuy/ρ) 0 ux uy

(uvx/ρ + vvy/ρ) 0 vx vy

. (5)

Se lo stato q e uguale alla soluzione stazionaria dell’Eq. (1), il termine a destra dell’ugualenell’Eq. (4) e nullo.Se il termine d′, che tiene conto dell’effetto viscoso dovuta alle variabili q′, viene trascurato,come usualmente fatto in acustica, si ottiene il sistema di equazioni di Eulero linearizzate

∂q′

∂t+ A

∂q′

∂x+ B

∂q′

∂y+ Cq′ = 0 . (6)

Per convenienza, l’Eq. (6) viene riscritta come

∂q′

∂t+

∂fx

∂x+

∂fy

∂y= s (7)

dove

fx =

u 0 ρ 00 u γp 00 1/ρ u 00 0 0 u

q′ fy =

v 0 0 ρ0 v 0 γp0 0 v 00 1/ρ 0 v

q′ (8)

rappresentano i flussi e il termine sorgente s

s =

0 0 0 00 (1− γ) (ux + vy) (γ − 1) px (γ − 1) py

− (uux + vuy) (1/ρ)x vy −uy

− (uvx + vvy) (1/ρ)y −vx ux

q′ (9)

rappresenta l’effetto dovuto al flusso medio non uniforme.

2 DISCRETIZZAZIONE SPAZIALE

La formulazione ai residui pesati dell’Eq. (7) puo essere scritta come∫

Ωφ

∂q

∂tdΩ +

∫

Ωφ · ∇~fdΩ =

∫

ΩφsdΩ (10)

dove con φ si indica una arbitraria funzione test e, per comodita, si usa la notazione q perindicare le grandezze perturbate q′.L’Eq. (10) viene discretizzata suddividendo il domino di calcolo Ω in un insieme di elementiE e limitando la scelta di q e φ a funzioni di tipo polinomiale definite all’interno di ognielemento. Non viene imposto alcun vincolo di continuita sulle funzioni q e φ, che risultano

A. Crivellini, F. Bassi 3

pertanto discontinue sulle interfacce degli elementi. Esprimendo l’integrale su Ω come sommadegli integrali sugli elementi E ed eseguendo una integrazione per parti, si ottiene la formulazionedebole dell’Eq (10). Se si considerano, per ogni elemento, funzioni φ che sono nulle all’esternodi E, si ottiene la equazione

∫

Eφ

∂q

∂tdΩ +

∫

∂Eφ~f · ~ndσ −

∫

E∇φ · ~fdΩ =

∫

EφsdΩ . (11)

A causa dell’approssimazione funzionale discontinua, la funzione di flusso nell’integrale di con-torno dell’Eq. (11) non e univocamente definita. Per assicurare un accoppiamento tra elementiadiacenti, che altrimenti mancherebbe completamente, il flusso fisico in direzione normale ~f · ~nviene sostituito da un flusso numerico h(ql, qr, ~n) che dipende dagli stati ql e qr sui due latidell’interfaccia e dal versore normale ~n. Gli stati ql e qr necessari per la valutazione del flussonumerico sono calcolati nei punti di integrazione sull’interfaccia utilizzando la rappresentazionepolinomiale della soluzione negli elementi adiacenti. Tra le varie funzioni di flusso numericodisponibili in letteratura, in questo lavoro e stato utilizzato il flusso di Roe

h (ql, qr, ~n) =12

(~f (ql) + ~f (qr)

)· ~n + P

[λ±

]P−1 (qr − ql)

, (12)

dove le matrici P e P−1 sono rispettivamente la matrice degli autovettori dello Jacobiano di~f (q) · ~n e la sua inversa e [λ±] e la matrice diagonale degli autovalori dello stesso. Tali matricisono calcolate in base allo stato medio di Roe, definito, indicando con H l’entalpia totale, dallerelazioni

ϑ =√

ρl√ρr +

√ρr

ρ =√

ρlρr; H = ϑHl + (1− ϑ) Hr (13)

u = ϑul + (1− ϑ) ur; v = ϑvl + (1− ϑ) vr .

Si noti che tale stato dipende esclusivamente dal valore delle grandezze del flusso medio sui duelati dell’interfaccia. Ovviamente, nel caso di flusso supersonico, h(ql, qr, ~n) coincide con il flussofisico in direzione normale calcolato in base allo stato sinistro (o destro).

3 INTEGRAZIONE TEMPORALE

L’integrazione temporale si basa su uno schema di Runge–Kutta esplicito (ERK) accurato alquarto ordine, proposto da A. Kennedy, H. Carpenter e Lewis, [12]. Tale schema e denomina-to Low Storage perche, a dispetto dell’ordine di accuratezza e del numero di stadi per passotemporale, richiede solo due registri per la memorizzazione dei dati. Questa caratteristica eovviamente apprezzabile in queste applicazioni che richiedono di memorizzare ingenti quantitadi dati. Applicando lo schema ERK all’equazione

M∂Q

∂t= R (t, q)

e considerando l’avanzamento della soluzione dall’istante t(n) all’istante t(n+1) si ha, per unoschma a s stadi,

4 56 Congresso Nazionale ATI

u0 = q(n)

ui = u0 + ai∆tM−1R(ti, ui

)+ ∆t

∑i−1

j=1bjM

−1R(tj , uj

), ti = t(n) + ci∆t

q(n+1) = u0 + ∆t∑s

j=1bjM

−1R(tj , uj

),

(14)

dove ai, bj e ci sono coefficienti propri di ogni schema, [12].

4 CONDIZIONI AL CONTORNO

L’imposizione delle condizioni al contorno viene realizzata costruendo, sulle facce di contorno,uno stato esterno che dipende dai dati imposti sul contorno e dalle informazioni uscenti daldominio di calcolo. Definito questo stato il calcolo del flusso viene effettuato attraverso il solutoredi Riemann utilizzato per gli elementi interni. Una conseguenza notevole e che, a differenza diquello che avviene usualmente negli schemi alle differenze finite ad elevato ordine di accuratezza,gli elementi di confine non hanno bisogno di alcun trattamento particolare.Inoltre, oltre agli ovvi benefici dovuti alla facilita d’implementazione e alla semplicita di gestionedel codice, non si introducono quelle modifiche che spesso, negli schemi classici, risultano essereuna delle maggiori fonti di errori e riduzione di accuratezza.

4.1 Condizioni al contorno di parete rigida

La condizione di parete rigida richiede che il flusso di massa che attraversa il contorno sia nullo.E possibile implementare tale condizione definendo uno stato esterno in base a una condizionedi simmetria locale

ρr = ρl

pr = pl

un,r = −un,l

ut,r = ut,l

(15)

e applicando poi il solutore di Riemann approssimato tra i due stati.

4.2 Condizioni al contorno non riflettenti

Uno dei requisiti principali in aeroacoustica computazionale e l’implementazione di condizionial contorno non riflettenti sulle frontiere di ingresso e uscita di flusso di massa (condizioni nonreflecting).Queste condizioni al contorno sono utilizzate per troncare un campo infinito in un dominio dicalcolo finito permettendo allo stesso tempo la propagazione delle perturbazioni ondose senzacausare indesiderate riflessioni non fisiche verso l’interno del dominio stesso.Diversi approcci sono presenti in letteratura tra i quali le condizioni non riflettenti di Giles, [8],le Radiation Boundary Conditions, [14], le tecniche di Perfectly Matched Layer (PML), [11], edi Buffer Layer (BL), [13].I primi due implicano la soluzione di un sistema di equazioni alle derivate parziali lungo ilcontorno mentre gli altri introducono attorno al dominio di calcolo vero e proprio una zona(buffer) in cui le equazioni sono modificate in maniera tale da annullare o attenuare le riflessionispurie verso l’interno del dominio (Figura 1).In questo lavoro viene utilizzata la tecnica di buffer layer e in tale zona il flusso medio acceleraprogressivamente fino a diventare supersonico sulle frontiere esterne del dominio allargato. Perottenere la voluta accelerazione del flusso medio basta modificare, nella zona di buffer layer,

A. Crivellini, F. Bassi 5

la componente di velocita normale alla frontiera sommando ad essa il valore ottenuto dallaespressione

ub = (uex − un)[ |dn (x, y)|

∆s

]m

se un < uex

ub = 0 se un ≥ uex

(16)

dove uex e il valore della velocita supersonica normale al contorno ed uscente da questo, ∆se lo spessore del buffer layer, dn(x, y) e la distanza dall’interfaccia tra buffer layer e campo dicalcolo del punto (x, y) ed m e un parametro che permette di controllare la legge di variazionedelle velocita. In base a precedenti esperimenti numerici, nell’Eq. (16) si e posto m = 2 e lavelocita uex corrisponde a un numero di Mach pari a 1.3. Si mette qui in evidenza che per poterapplicare tale tecnica il flusso deve diventare uscente in corrispondenza dei contorni esterni deldominio allargato; nel buffer layer adiacente a una frontiera di ingresso si ha quindi l’inversionedella direzione della velocita del flusso medio.Nella zona formata dall’intersezione di due buffer layer adiacenti a contorni che formano unospigolo tra di loro, l’Eq. (16) viene utilizzata una volta per ciascuno di essi, ottenendo il risultatomostrato in Figura 2 (flusso medio interno in direzione x con numero di Mach uguale a 0.5).

x

y

-150 -100 -50 0 50 100 150-150

-100

-50

0

50

100

150

dominio interno

buffer layer

Figura 1: Esempio di buffer layer

x

y

80 90 100 110 120 130

80

100

120 buffer layer

dominio interno

Figura 2: Buffer layer in presenza di spigoli

E importante mettere in evidenza gli effetti che la spaziatura della griglia di calcolo nella zona dibuffer layer puo avere sul passo di integrazione temporale. Se la dimensione caratteristica deglielementi, dr, fosse uniforme in tutto il dominio di calcolo, gli elementi piu vicini alla frontieraesterna della zona di buffer layer porrebbero un limite piu stringente al valore massimo del passodi integrazione temporale: infatti, ∆t = CFL

[dr

/(c +

∣∣~v∣∣)] risulterebbe minimo dove c +∣∣~v∣∣ e

massimo. Per evitare il conseguente aumento di tempo di calcolo complessivo, e opportuno chela griglia di calcolo presenti elementi di dimensioni maggiori nella zona di buffer layer.

5 ESEMPIO DI APPLICAZIONE

Il codice di calcolo sviluppato e stato provato con il problema della categoria 5 del “ThirdComputational Aeroacoustics (CAA) Workshop on Benchmark Problems”. Questo test case sipropone di modellare la presenza di una sorgente acustica all’interno di un getto bidimensionale e

6 56 Congresso Nazionale ATI

di verificare la capacita del metodo di calcolare l’esatta propagazione acustica in un flusso medionon uniforme caratterizzato dalla presenza di un forte gradiente. Lo stato medio e rappresentatodalle relazioni

u (y) =

uj se 0 ≤ y < h

u∞ + (uj − u∞)− ln 2( yb−h

b )2

se y ≥ h

1ρ (y)

= −12

γ − 1yp

[u (y)− uj ] [u (y)− u∞] +1

ρ∞[u (y)− u∞](uj − u∞)

+1

ρ∞[uj − u (y)](uj − u∞)

Il valore della pressione statica e invece costante e il getto e considerato simmetrico e paralleloall’asse delle x. Per descrivere la sorgente acustica, al secondo membro dell’Eq. (7) viene aggiuntoil termine

r = (0, r2, 0, 0)T ; r2 = A−B ln 2(x2+y2) cosωt. (17)

Le Tabelle 1, 2 e 3 riportano i valori dimensionali delle grandezze che caratterizzano il problema.

Tabella 1: Stato termodinamico

p = 101330 [Pa] p∞ = 101330 [Pa]Tj = 166.6667 [K] T∞ = 300 [K]uj = 517.4569 [m/s] u∞ = 0 [m/s]

Tabella 2: Costanti del gas ideale

γ = 1.4R = 286.8875 [m2/(s2K)]

Tabella 3: Costanti geometriche

l = h + b = 1.0 [m]h = 0.6 [m]b = 0.4 [m]

Tabella 4: Costanti del termine sorgente

A B St ω T[Kg/(ms3)] [1/m2] [rad/s] [s]

caso A 0.001 8 0.14 227.589 2.761× 10−2

caso B 0.001 8 0.60 1950.767 6.442× 10−3

Tali valori determinano un numero di Mach nella parte centrale del getto pari a 2.Come si vede dalla Tabella 4 il problema richiede simulazioni per due diversi numeri di Strouhal(St = fL/uj), calcolati in base alla velocita del getto e alla distanza L di mezza velocita1, cuiovviamente corrispondono due diversi valori di ω (ω = 2πStuj/L) nell’Eq. (17).

1L = distanza tra i punti per i quali la velocita e la meta della massima del getto

A. Crivellini, F. Bassi 7

Il dominio fisico prescritto e di −5 m ≤ x ≤ 50 m e 0 ≤ y ≤ 10 m e si richiede di calcolare ilsuono radiato, come p2, lungo le linee a y = 10 m e x = 50 m, oltre al valore della pressioneacustica, all’istante iniziale di un ciclo, per y = 1 m.

5.1 Procedura e griglie di calcolo

Tutte le grandezze sono state adimensionalizzate in base ai valori della densita e della velocita delsuono dello stato indisturbato e della lunghezza l = L/2 = 1 m. Si noti che questo determinaun valore della ampiezza della sorgente acustica dell’Eq. (17) pari a A ' 2.030736 × 10−11.Questo aspetto rende ancora piu evidente la necessita di separare esplicitamente le variabilifluidodinamiche da quelle acustiche. Come si vedra piu avanti, infatti, si ottengono valori delleperturbazioni che se sommate allo stato mediato sarebbero confrontabili se non inferiori aglierrori di troncamento dovuti alla precisione del calcolatore.In tutte le simulazioni la soluzione avanza nel tempo con un valore di CFL massimo pari a 0.6e si basa su una discretizzazione spaziale con polinomi di terzo grado (P3) cui corrisponde unaaccuratezza spaziale del quarto ordine.Sono state utilizzate quattro diverse griglie di calcolo formate da elementi triangolari e dotatedi un buffer layer di spessore (∆s/l) pari a 20. In tutti i casi il dominio di calcolo e suddivisoin zone con elementi uguali all’interno di ogni zona. I parametri fondamentali di tali grigliesono riportati nalla Tabella 5 dove sono indicati il numero di elementi compresi in ogni trattoe il numero di elementi totali. La Figura 3 mostra, a titolo di esempio, una griglia di calcolomolto piu rada di quelle effetivamente utilizzate. In particolare si puo notare come gli elementisiano addensati in prossimita del gradiente di flusso medio (0.6 < y/l ≤ 1.7), diradati all’internodel buffer layer e come i passaggi tra le varie zone siano bruschi senza alcuna transizione. Glielementi delle griglie hanno forma di triangolo rettangolo, con un valore massimo del rapportotra le dimensioni dei cateti pari a 14.

Tabella 5: Griglie di calcolo

griglia −25 ≤ x < −5 −5 ≤ x < 55 55 ≤ x < 70 0 ≤ y < 0.6 0.6 ≤ y < 1.7 1.7 ≤ y < 15 15 ≤ y < 30 Nelem

G1 10 70 10 3 18 13 10 7920G2 13 91 13 5 24 17 13 13806G3 15 105 15 7 27 20 15 18630G4 20 140 20 7 27 26 20 28800

x

y

-25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 60 65 700

5

10

15

20

25

30

Figura 3: Tipica griglia di calcolo

8 56 Congresso Nazionale ATI

5.2 Caso A (St = 0.14)

Il primo caso e quello caratterizzato dalla minore frequenza della sorgente acustica cui cor-risponde una maggiore lunghezza d’onda delle perturbazioni. Per questo motivo i risultati ot-tenuti con la griglia piu rada (G1) appaiono di buona qualita e, come si nota nelle Figure 4, 5 e 6,sono praticamente coincidenti con quelli calcolati sulle griglie piu fini. Solo nel tratto compresotra 0 ≤ y/l < 1 (parte centrale del getto) si puo notare una lieve differenza tra i risultati ot-tenuti nelle varie simulazioni. Come si e gia osservato, le Figure mettono in evidenza i valoriparticolarmente piccoli delle perturbazioni acustiche.Per permettere che il dato iniziale del calcolo non lasci traccia nella soluzione calcolata, lasimulazione e stata fatta avanzare per 15 periodi della sorgente corrispondenti ad un tempofisico di circa 0.414 s. Si e riscontrata una non perfetta periodicita del suono radiato (p2) soloin una determinta zona della sezione a y/l = 10 (in corrispondenza dei valori piu piccoli delleperturbazioni acustiche). La Figura 7 mostra i valori calcolati per tale zona in alcuni periodisuccessivi.Infine, la Tabella 6 riassume le risorse di calcolo richieste usando un solo processore CompaqAlpha ev67 a 667 MHz.

x

p2

0 10 20 30 40 50

10-14

10-13

10-12

G1G2G3

Figura 4: Suono radiato a y = 10 (St =0.14)

y

p2

0 2 4 6 8 10

10-12

10-11

10-10

G1G2G3

Figura 5: Suono radiato a x = 50 (St =0.14)

Tabella 6: Risorse di calcolo (St = 0.14)

griglia memoria richiesta numero di cicli cpu time[Mbyte] [s]

G1 50 13180 2.52× 104

G2 87 17676 5.82× 104

G3 134 19940 8.91× 104

5.3 Caso B (St = 0.60)

Essendo la lunghezza d’onda caratteristica di questo problema piu piccola della precedente sonostate utilizzate le tre griglie di calcolo con maggiore risoluzione. Il tempo totale di simulazione,

A. Crivellini, F. Bassi 9

x

p

0 10 20 30 40 50

-6.0x10-06

-3.0x10-06

0.0x10+00

3.0x10-06

6.0x10-06G1G2G3

Figura 6: Perturbazione acustica a y=1all’inzio di un ciclo (St = 0.14)

xp2

10 20

1.0x10-14

2.0x10-14

3.0x10-14

4.0x10-14

5.0x10-14

6.0x10-147.0x10-14

T=15T=16T=17T=18T=19

Figura 7: Suono radiato a y=10 in periodisuccessivi (St = 0.14)

x

y

-20 -10 0 10 20 30 40 50 60 700

5

10

15

20

25

30

p: -7.10x10-07 -5.07x10-07 -3.04x10-07 -1.01x10-07 1.01x10-07 3.04x10-07 5.07x10-07 7.10x10-07

Figura 8: Perturbazione acustica all’inzio di un ciclo (St = 0.14)

10 56 Congresso Nazionale ATI

corrispondente a 60 periodi della sorgente, e di 0.387 s.Lievi differenze tra le soluzioni sono percepibili solo nella sezione x/l = 50 (Figura 10) situataa maggior distanza dalla sorgente e in prossimita della zona di buffer layer. La Figura 9 mostrainvece che le soluzioni a y/l = 10 sono praticamente coincidenti tranne che, per la griglia G2,nella zona adiacente al confine del dominio fisico.Per quanto riguarda la pressione acustica valutata all’inizio del periodo, a y/l = 1, si puonotare anche in questo caso che la griglia di minore costo computazionale ha una soluzioneassai prossima alle altre due (che si ricoprono in maniera pressoche perfetta), con un andamentogenerale della soluzione del tutto analogo, ma con picchi leggermente smussati a causa dellamaggiore dissipazione numerica.La Tabella 7, infine, riporta le risorse di calcolo impiegate per le simulazioni.

x

p2

0 10 20 30 40 50

10-16

10-15

10-14

G2G3G4

Figura 9: Suono radiato a y = 10 (St =0.60)

y

p2

0 2 4 6 8 1010-16

10-15

10-14

10-13

G2G3G4

Figura 10: Suono radiato a x = 50 (St =0.60)

Tabella 7: Risorse di calcolo (St = 0.60)

griglia memoria richiesta numero di cicli cpu time[Mbyte] [s]

G2 87 16497 5.44× 104

G3 134 18611 8.22× 104

G4 206 18705 1.30× 105

CONCLUSIONI

In questo lavoro e stato sviluppato un metodo agli elementi finiti discontinui di Galerkin perlo studio della propagazione acustica descritta dalle equazioni di Eulero linearizzate attorno aduno stato medio non uniforme.La discretizzazione spaziale agli elementi finiti discontinui consente di ottenere elevata accuratez-za spaziale anche su griglie di calcolo fortemente discontinue. Lo schema di Runge–Kutta esplic-ito a cinque passi permette di avanzare la soluzione nel tempo con elevato ordine di accuratezzae richiede solo due registri per la memorizzazione.

A. Crivellini, F. Bassi 11

x

p

10 20 30 40 50

-3.0x10-07

-2.0x10-07

-1.0x10-07

0.0x10+00

1.0x10-07

2.0x10-07

G2G3G4

Figura 11: Perturbazione acustica a y = 1 all’inzio di un ciclo (St = 0.60)

x

y

-20 -10 0 10 20 30 40 50 60 700

5

10

15

20

25

30

p: -1.40x10-07 -7.78x10-08 -1.56x10-08 4.67x10-08 1.09x10-07

Figura 12: Perturbazione acustica all’inzio di un ciclo (St = 0.60)

12 56 Congresso Nazionale ATI

Per trattare le frontiere attraversate da flusso di massa e stata utilizzata una tecnica di bufferlayer. Tale tecnica e risultata di facile implementazione e particolarmente efficace per ridurre leriflessioni spurie dalle frontiere esterne del dominio di calcolo.Gli esempi di calcolo riportati confermano le notevoli caratteristiche di accuratezza del metodo.

Riferimenti bibliografici

[1] H. Atkins. Application of the Discontinuous Galerkin method to Acoustic Scatter Problems.In C. K. W. Tam and J. C Hardin, editors, Second Computational Aeroacoustics (CAA)Workshop on Benchmark Problems, volume NASA Conference Publication 3352, pages1300–1400, June 1997.

[2] H. Atkins. Continued Development of Discontinuous Galerkin Method for ComputationalAeroacoustic Application. AIAA Paper, 97-1581, Maggio 1997.

[3] H. Atkins and D. Lockard. A High-Order Method Using Unstructured Grids for theAeroacoustic Analysis of Realistic Aircraft Configuration. AIAA Paper, 99-1945, 1999.

[4] H. Atkins and D. Lockard. An Application of the Quadrature Free Discontinuous GalerkinMethod. xxxx, 1999.

[5] H. Atkins and C. Shu. Quadrature-Free Implementation of Discontinuous Galerkin Methodfor Hyperbolic Equation. AIAA Journal, 36(5):775–782, 1997.

[6] B. Cockburn, S. Hou, and C. W. Shu. TVB Runge-Kutta Local Projection Discontin-uous Galerkin Finite Element Method for Conservation Laws IV: The MultiDimensionalEstension. Mathematics of Computation, 54(190):545–581, 1990.

[7] B. Cockburn and C. W. Shu. TVB Runge-Kutta Local Projection Discontinuous GalerkinFinite Element Method for Conservation Laws II: General Framework. Mathematics ofComputation, 52(186):411–435, 1989.

[8] M. Giles. Non Reflecting Boundary Condition For Euler Equation Calculation. AIAAJournal, 28(12):2050–2058, 1990.

[9] J. C Hardin, J. R. Ristorcelli, and C. K. W. Tam. ICASE/Larc Workshop on BenchmarkProblems in Computational Aeroacoustics (CAA). NASA Conference Publication 3300,May 1995.

[10] R. Hixon, S. Shih, and R. Mankbadi. Evalutation of Boundary Condition for ComputationalAeroacoustics. AIAA Journal, 33(11):2006–2012, 1995.

[11] F. Hu. On Absorbing Boundary Condition for Linearized Euler Equation by a PerfectlyMatched Layer. Journal of Computational Physics, 129:201–219, 1996.

[12] A. Kennedy, H. Carpenter, and R. Lewis. Low-Storage, Explicit Runge-Kutta Schemes forthe Compressible Navier Stokes Equation. ICASE Report, NASA/CR-1999-209349(99-22),1999.

[13] D. M. Nark. The Use of Staggered Schemes and an Absorbing Buffer Zone for Computation-al Aeroacoustics. In J. C Hardin, J. R. Ristorcelli, and C. K. W. Tam, editors, ICASE/LarcWorkshop on Benchmark Problems in Computational Aeroacoustics (CAA), volume NASAConference Publication 3300, pages 233–244, May 1995.

[14] C. K. W. Tam. Computational Aeroacoustics: Issues and Methods. AIAA Journal,33(10):1788–1795, 1995.

A. Crivellini, F. Bassi 13

[15] C. K. W. Tam and J. C Hardin. Second Computational Aeroacoustics (CAA) Workshopon Benchmark Problems. NASA Conference Publication 3352, June 1997.

[16] C. K. W. Tam and J. C. Web. Dispersion-Relation-Preserving Schemes for ComputationalAcoustics. journal of Computational Physics, 107:262–281, 1993.