J Diaz

P Joly

A time domain analysis of PML models in acoustics

Keywords: Perfectly Matched Layers, Acoustics, Stability, Error Estimates, Cagniard-de Hoop. PACS:

We present in this paper a time-domain analysis of PML's for non-advective and advective acoustics. We focus our attention on time-stability and error estimates (with respect to the parameters of the layers). The main new technical tool is the Cagniard-de Hoop method. Our theoretical results are validated and illustrated by various numerical results.

In section 1, we first derive an error analysis for the PML for non advective acoustics. This work is a continuation of [START_REF] Diaz | An analysis of higher boundary conditions for the wave equation[END_REF] that concerned absorbing boundary conditions. The main point is the derivation of an explicit solution using the Cagniard-De Hoop technique [START_REF] Cagniard | Reflection and refraction of progressive seismic waves[END_REF][START_REF] De Hoop | The surface line source problem[END_REF][START_REF] De Hoop | Analytic time-domain performance analysis of absorbing boundary conditions and perfectly matched layers[END_REF]. Our results show the exponential accuracy of PML's with respect to the damping factor σ and the width L of the layer. Section 2 is devoted to classical PML's for advective acoustics: the propagation of waves in the presence of a uniform mean flow is modeled by linearized Euler equations. We revisit the well-known instability of standard PML's by the Cagniard-De Hoop technique.

In section 3, we are concerned by the question of the stabilization of PML's for advective acoustics when the flow is orthogonal to the artificial boundary. We proposed a stable PML (already introduced in [START_REF] Diaz | Stabilized perfectly matched layer for advective acoustics[END_REF]) which consists, in a first step, of sticking to the physical domain a non advective propagation model designed in such a way that no reflection is produced at the interface between the physical domain and the artificial one (this is the "PML property"), and in a second step, of using standard PML's for the artificial domain. We derive the stability and accuracy analysis for this new PML.

Our theoretical results are illustrated by various numerical results. Of course the present work has to be compared with some recent ones in the literature, in particular by Hagström and Nazarov [START_REF] Hagström | Absorbing layers and radiation boundary condition for jet flow simulations[END_REF], Abarbanel, Gottlieb and Hesthaven [START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF] and more specially by Hu [START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF] whose work on stabilized PML for advective acoustics appears to be close to what we are doing in this paper as it will be briefly discussed in the last section of this paper.

An appendix is devoted to a partial solution to the problem of constructing stable PML's in the case of an oblique (with respect to the boundaries of the rectangular computational domain) mean flow.

PML for acoustics

Let us consider the two dimensional wave equation

1 c 2 ∂ 2 u ∂t 2 -∆u = f (x, t), x = (x 1 , x 2 ) ∈ IR 2 , t > 0.
(

) 1
If the initial data of the problem and the domain of interest are supported in the left half-space IR 2 -= IR -× IR it is natural to try to reduce the computation to this half-space. There are two main classes of method to do so. The first ones are the so-called absorbing boundary condition, whose principle consists in imposing a specific condition on the boundary designed so that the waves leave the computational domain without being (too much) reflected. Another way to bound the computational domain consists in adding a vertical layer where [START_REF] Bérenger | A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF] is modified by an absorption term (see figure 1). In [START_REF] Bérenger | A Perfectly Matched Layer for the Absorption of Electromagnetic Waves[END_REF][START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF], Bérenger proposed to introduce a very special anisotropic damping term acting only in the direction which is orthogonal to the interface. An alternative to the original split form of its model (see section 2) consists in replacing the derivation with respect to x 1 by the operator D σ

x 1 formally defined by :

D σ x 1 = ∂ ∂t + σ(x 1 ) -1 ∂ ∂t ∂ ∂x 1 ,
where σ is a non-negative function of x 1 .

More rigorously ψ = D σ x 1 φ means :

∂ ∂t + σ(x 1 ) ψ = ∂ ∂t ∂ ∂x 1 φ and ψ(t = 0) = 0.
The equation in the PML then becomes :

1 c 2 ∂ 2 u ∂t 2 -(D σ x 1 ) 2 u - ∂ 2 u ∂x 2 2 = 0, x = (x 1 , x 2 ) ∈ [0 ; L] × IR, t > 0. ( 2
)
The most interesting property of the PML is that a wave propagating in the computational domain is transmitted to the absorbing layer without being reflected. As this method is moreover really easy to implement (even in the corners of the domain), it has rapidly attracted a lot of people in different fields of application.

The Main results

The case of an infinite PML. We are first interested in the explicit expression of the fundamental solution u σ,∞ of the wave equation in the left half space IR 2 -with an infinite PML in the right half-space. u σ,∞ satisfies :

          
Find u : IR 2 × IR + → IR, null for t < 0 such that

1 c 2 ∂ 2 u ∂t 2 -(D σ x 1 ) 2 u - ∂ 2 u ∂x 2 2 = δ(x -x S ) × δ(t), in IR 2 × IR + . (3)
where x S = (-h, 0), σ(x 1 ) = 0 if x 1 ≤ 0 and σ(x 1 ) > 0 if x 1 > 0.

Let us define r(x) = |xx S | and the function θ(x) by :

θ(x) ∈ ]0, 2π] , x -x S = ( r(x) cos θ(x), r(x) sin θ(x) ) t ,
and the function Σ(x 1 ) by : Σ(x 1 ) =

x 1 0 σ(x) dx. The functions A(x, t) and B(x, t) are defined by (r = r(x) and θ = θ(x)) :

A(x, t) = | cos θ| Σ(x 1 ) t r > 0 and B(x, t) = | sin θ| Σ(x 1 ) t 2 r 2 - 1 c 2 .
Theorem 1 The expression of the solution

u σ,∞ (x, t) = G σ,∞ i (x, t) of prob- lem (3) is given by : G σ,∞ i (x, t) = H(c t -r(x)) 2π t 2 -r(x) 2 c 2 , x 1 ≤ 0 G σ,∞ i (x, t) = H(c t -r(x)) 2π t 2 -r(x) 2 c 2 e -A(x,t) cos [B(x, t)] , x 1 > 0. ( 4
)
where H denotes the Heaviside function

Remark 2

The restriction of G σ,∞ i (x, t) to the left half-space IR 2 -does not depend on σ and is nothing more than the restriction to the left half-space of the fundamental solution of the wave equation in the whole plane : the infinite PML does not produce any parasitical reflection.

The case of a finite PML. We now consider a PML of finite width L with a Neumann condition at its bottom end. Problem (3) then becomes :

                       Find u : IR 2 × IR + → IR, null for t < 0 such that 1 c 2 ∂ 2 u ∂t 2 -(D σ x 1 ) 2 u - ∂ 2 u ∂x 2 2 = δ(x -x S ) × δ(t), in [-∞ ; L] × IR × IR + , ∂u ∂x 1 (L, x 2 , t) = 0.
(5) We define the image point source, symmetrical to the point source with respect to the line x 1 = L by :

x * S = (h + 2L, 0), we denote r * (x) = |xx * S | and we define the function θ * (x) by :

θ * (x) ∈ ] π 2 , 3π 2 ] , x -x * S = ( r * (x) cos θ * (x), r * (x) sin θ * (x) ) t .
x * S = (h + 2L, 0)

r ( x ) r * ( x ) θ * (x, y) x 1 x 2 θ x S = (-h, 0) x 1 = L Fig. 2

. Illustration of the notation

We also define the function Σ * (x 1 ) by :

Σ * (x 1 ) = L 0 σ(x) dx + L x 1 σ(x) dx (= 2Σ(L) for x 1 < 0),
and the functions A * (x, t) and B * (x, t) by ( r * = r * (x) and θ * = θ * (x)) :

A * (x, t) = | cos θ * |Σ * (x 1 ) t r * > 0 and B * (x, t) = | sin θ * |Σ * (x 1 ) t 2 r * 2 - 1 c 2 .
Theorem 3 The expression of the solution u σ,L (x, t) = G σ,L (x, t) of problem (5) is given by :

G σ,L (x, t) = G σ,∞ i (x, t) + G σ,L r (x, t) with G σ,L r (x, t) = H(c t -r * (x)) 2π t 2 -r * (x) 2 c 2 e -A * (x,t) cos [B * (x, t)] . (6)
This theorem can easily be deduced from theorem 1 thanks to the image principle. Therefore, only theorem 1 will be proved in this paper.

Error Estimates. Let us consider the approximation, in the upper halfspace, of the solution u of problem (1) (the wave equation with a "regular" point source), by the function u σ,L , solution of problem :

                       Find u σ,L : IR 2 × IR + → IR, null for t < 0 such that 1 c 2 ∂ 2 u σ,L ∂t 2 -(D σ x 1 ) 2 u σ,L - ∂ 2 u σ,L ∂x 2 2 = δ(x -x S )×f (t), in [-∞; L] × IR × IR + , ∂u σ,L ∂x 1 (L, x 2 , t) = 0.
(7) After having defined the function Φ(t), null for t < (h + 2L)/c, by :

Φ(t) = Log   t + t 2 -((2L + h)/c) 2 ((2L + h)/c)   , if h + 2L c ≤ t ≤ h + 2L c + T, Φ(t) = Log   t + t 2 -(t -T ) 2 t -T   , if t > h + 2L c + T.
we can state the Theorem 4 One has the uniform estimates :

(u -u σ,L )(., t) L ∞ (IR 2 -) ≤ e -2 2L+h c 2 t Lσ 2π Φ(t) f L ∞ . ( 8
)
where σ = 1 L L 0 σ(x)dx is the mean value of σ in the PML.

This result appeals the following comments :

• The error converges spectrally to 0 (in the uniform norm) when σ (or L) goes to infinity. • For given t, the upper bound in the estimate diminishes when the distance h from the source to the absorbing layer increases. This is coherent with the physical intuition and numerical observations. • Concerning the behavior of the error for large t, if we assume that T < +∞, we observe that the right hand side in the estimate behaves for large t as

1 2π 2T t ,
which shows that, for all σ, L and h, the error converges uniformly to 0 when t tends to +∞.

• As the plane wave analysis might suggest, the shape of the function σ(x 1 ) does not influence the estimates, only its mean value between 0 and L has an importance. • We could have imagined, once again from the plane wave analysis, that σ and L did play the same role (i.e. increasing L is equivalent to increasing σ), but the uniform estimates show that increasing L allows us, not only to increase the absorption thanks to the term L σ, but also to move away the point source from the absorbing layer thanks to the term 2L + h.

The case of a rectangular domain surrounded by PML. One of the advantages of PML when compared the to absorbing boundary conditions (ABC's) is the simplicity of their implementation in the corner. By the principle of images and Cagniard-de Hoop method, it is also easy to compute the analytical solution in the case of a rectangular domain surrounded by PML's (whereas it is not possible with ABC's).

We now suppose that the computational domain Ω is the square [-h ; h] 2 and is surrounded by PML's of width L. We impose a Neumann condition at the end of the layers. We then consider the problem :

                                             Find u : IR 2 × IR + → IR such that 1 c 2 ∂ 2 u ∂t 2 -(D σ x 1 ) 2 u -(D σ x 2 ) 2 u = δ(x) × δ(t), in [-L; L] 2 × IR + , ∂u ∂x 1 (x, t) = 0, for x 1 = ±L ∂u ∂x 2 (x, t) = 0, for x 2 = ±L, u(x, t) = 0, for t < 0. ( 9
)
Let us define the image points source (x ij ) (i,j)∈Z 2 as

x ij = (2i(h + L), 2j(h + L)).
We set r ij (x) = |xx ij | and we define the functions θ ij (x) by :

θ ij (x) ∈ ] -π, π] , x = ( r ij (x) sin θ ij (x), r ij (x) cos θ ij (x) ) t ,
and the non-negative functions (Σ i (x)) i∈Z by :

                     Σ i (x) = (2i -1) L 0 σ(x) dx + x -L σ(x) dx if i < 0, Σ i (x) = (2i -1) L 0 σ(x) dx + L x σ(x) dx if i > 0, Σ 0 (x) = x 0 σ(x) dx .
Finally we define the functions A ij (x, t) and B ij (x, t) by :

A ij (x, t) = (| cos θ ij (x)|Σ j (x 2 ) + | sin θ ij (x)|Σ i (x 1 )) t r ij (x)
> 0 and

B ij (x, t) = (|cos θ ij (x)| Σ i (x 1 ) + |sin θ ij (x)| Σ j (x 2 )) t 2 r ij (x) 2 - 1 c 2 .
x 00

x 01 x 10

x 11

x 1-1

x -10

x i0

x i1

x i-1

x ij

x i-j

x -i0

x -i1

x -i-1

x -ij

x -i-j

x 0j

x 0-j

x -1j

x -1-j

x 1j

x 1-j

x -11

x -1-1

2(L + h) 2(L + h) h L x 0-1 Fig. 3. The image points source Theorem 1.1 The solution u σ,h,L (x, t) = G σ,h,L (x, t) of problem (9) is given by : G σ,h,L (x, t) = (i,j)∈Z 2 G σ,h,L i,j (x, t)
where

G σ,h,L ij (x, t) = H(c t -r ij (x)) 2π t 2 - r ij (x) 2 c 2 e -A ij (x,t) cos [B ij (x, t)] . (10)
Remark 5 For given t and x, G σ,h,L (x, t) is a finite sum. Indeed if i and j are such that r ij (x) < ct, then G σ,h,L ij (x, t) is null.

We do not present here the details of the proof which is very similar to the half-plane case, we refer the reader to [START_REF] Diaz | Approches analytiques et numériques de problèmes de transmission en propagation d'ondes en régime transitoire. application au couplage fluidestructure et aux méthodes de couches parfaitement adaptées[END_REF] for details.

Proof of the theorems

Expression of the Green's function. Let G σ,∞ i (x, t) be the solution of problem [START_REF] Zhao | A General Approach for the Deveopment of Unsplit-Field Time-Domain Implementations of Perfectly Matched Layers for FDTD Grid Truncation[END_REF]. The first step of the Cagniard-de Hoop's method consists in applying successively to G σ,∞ i (x, t) the Laplace transform in time and the Fourier transform in the space variable x 2 (we denote by s (resp. k) the dual variable of t (resp. x 2 )). The resulting function x 1 → G σ,∞ i (x 1 , k, s) satisfies the ordinary differential equation :

- s s + σ(x 1 ) d dx 1 ( s s + σ(x 1 ) d G σ,∞ i dx 1 ) + (k 2 + s 2 c 2 ) G σ,∞ i = δ(x 1 + h), Introducing G X 1 (x 1 , s) = x 1 + 1 s x 1 0 σ(x)dx, k, s = G σ,∞ i (x 1 , k, s) we get : - d 2 G dX 2 1 + (k 2 + s 2 c 2 ) G = δ(X 1 + h).
which leads to the expressions :

G(X 1 , k, s) = e -(k 2 + s 2 c 2 ) 1 2 |X 1 +h| 2(k 2 + s 2 c 2 ) 1 2 =⇒ G σ,∞ i (x 1 , k, s) = e -(k 2 + s 2 c 2 ) 1 2 x 1 +h+ Σ(x 1 ) s 2(k 2 + s 2 c 2 ) 1 2
.

Let us remark that, if x 1 is negative, Σ(x 1 ) = 0 and

G σ,∞ i (x 1 , k, s) = e -(k 2 + s 2 c 2 ) 1 2 |x 1 +h| 2(k 2 + s 2 c 2 ) 1 2
, which is nothing but the Fourier-Laplace transform of the fundamental solution of the wave equation :

G σ,∞ i (x 1 , x 2 , t) = H(c t -r(x)) 2π t 2 -r(x) 2 c 2 , x 1 < 0.
Let us now consider x 1 > 0 so that

x 1 + h + Σ(x 1 ) s > 0 and G σ,∞ i (x 1 , k, s) = e -(k 2 + s 2 c 2 ) 1 2 x 1 +h+ Σ(x 1 ) s 2(k 2 + s 2 c 2 ) 1 2
.

Applying the x 2 inverse Fourier transform to G σ,∞ i we get :

G σ,∞ i (x 1 , x 2 , s) = 1 4π +∞ -∞ e -(k 2 + s 2 c 2 ) 1 2 x 1 +h+ Σ(x 1 ) s -ikx 2 (k 2 + s 2 c 2 ) 1 2 dk.
Contrary to what a reader familiar with the Cagniard-de Hoop could expect, the term in the exponential is not an homogeneous function in (k, s). Applying the usual change of variable k = ps/c allows us however to obtain :

G σ,∞ i (x 1 , x 2 , s) = e -Σ(x 1 ) c (1+p 2 ) 1 2 e -s c [(1+p 2 ) 1 2 (x 1 +h)+ipx 2 ] 4π(1 + p 2 ) 1 2 dp, ≡ +∞ -∞ Ψ(p) dp .
We then have to deal with the product of two exponential functions. The first one does not depend on s and the other one is the one that usually appears in most applications of the Cagniard-de Hoop method.

Let us introduce the so-called Cagniard-de Hoop contour Γ defined as :

Γ = Γ + ∪ Γ -with Γ ± =    p = γ ± (t) ≡ -i ct r sin θ ± cos θ c 2 t 2 r 2 -1, t ≥ r c    . so that (1 + p 2 ) 1 2 (x 1 + h) + ipx 2 = ct ∈ IR + for all p ∈ Γ.
It is clear that the two curves Γ ± are symmetric the one from the other with respect to the imaginary axis and meet at point -i sin θ (for t = r/c). Actually Γ is nothing but the branch of the hyperbola of equation:

Y 2 sin 2 θ - X 2 cos 2 θ = 1, (p = X + iY, (X, Y ) ∈ IR 2 ),
which is located in the upper half-space Y = ℑm p > 0 if sin θ < 0. Note that this hyperbola does not intersect the two branch cuts of Ψ. All this information is summarized in figure 4.

Let us denote by D the real line and by Ω the connected part of the complex plane delimited by D and Γ. Let ρ > 0 a parameter devoted to tend to +∞. We set:

D ρ = {p ∈ D / |p| ≤ ρ} , Γ ρ = {p ∈ Γ / |p| ≤ ρ} and C ρ = {p ∈ Ω / |p| = ρ} . Note that D ρ ∪ C ρ ∪ Γ ρ is a closed curve (see figure 5). Since Ψ(p) is analytic D ℑm(p) ℜe(p) Γ - Ω Γ + i -i sin θ Fig. 4. The contours Γ and D -ρ ρ D C -(ρ) C + (ρ) Γ -(ρ) Γ + (ρ) i Fig. 5. The closed contour D ρ ∪C ρ ∪Γ ρ
in Ω, the integral of Ψ along D ρ ∪ C ρ ∪ Γ ρ , that we orient according to figure 5, is identically 0 :

Dρ Ψ(p) dp + Cρ Ψ(p) dp + Γρ Ψ(p) dp = 0.
Thanks to the choice of the square root and since x 1 + h > 0, the function Ψ(p) decays exponentially to 0 when ℑm p goes to +∞. As a consequence, it is easy to show that (Jordan's lemma) :

lim ρ→+∞ Cρ Ψ(p) dp = 0, which implies +∞ -∞ Ψ(p) dp = - Γ Ψ(p) dp and G σ,∞ i (x 1 , x 2 , s) = - 1 4π Γ e -Σ(x 1 ) c (1+p 2 ) 1 2 e -s [ (1+p 2 ) 1 2 x 1 +h c +ip x 2 c ] (1 + p 2 ) 1 2
dp.

We use the parametrisations p = γ + (t) and p = γ -(t), for t ≥ r/c, respectively along Γ + and Γ -and remark that :

• (1 + p 2 ) 1 2 ( x 1 + h c ) + ip x 2 c = t, ( by
construction),

• dp

(1 + p 2 ) 1 2 = ± dt (t 2 -r 2 c 2 ) 1 2
, on Γ ± ,

• (1 + p 2 ) 1 2 c = t r 2 (x 1 + h) ∓ i x 2 r 2 (t 2 - r 2 c 2 ) 1 2 , on Γ ± .
Therefore, since t goes from +∞ to r c on Γ + and from r c to +∞ on Γ -:

G σ,∞ i (x 1 , x 2 , s) = 1 2π +∞ r c e -(x 1 +h)Σ(x 1 ) t r 2 (t 2 -r 2 c 2 ) 1 2 cos x 2 r 2 Σ(x 1 )(t 2 - r 2 c 2 )
1 2

e -st dt.

We can now conclude, using the injectivity of the Laplace-Fourier transform.

Error Estimates Let u and u σ,L be the respective solutions of ( 1) and [START_REF] Fauqueux | Eléments finis mixtes spectraux et couches absorbantes parfaitement adaptées pour la propagation d' ondes élastiques en régime transitoire[END_REF]. We introduce the error (or reflected field) defined as e σ,L = u σ,L -u. Using the fact that f is supported in [0, T ] and

G σ,L r in [r * /c, ∞],
we obviously have e σ,L (x, t) = 0 if t ≤ r * /c and :

e σ,L (x, t) = +∞ max(r * /c, t-T ) G σ,L r (x, τ ) f (t -τ ) dτ if t > r * c . (11)
We deduce that :

           | e σ,L (x, t)| ≤ f L ∞ (0,t) • G σ,L r (x, .) L 1 ( r * c ,t) if r * c ≤ t ≤ r * c + T, | e σ,L (x, t)| ≤ f L ∞ (0,T ) • G σ,L r (x, .) L 1 (t-T,t) if t > r * c + T. (12)
Let us now introduce the two disjoint sets Ω 1 (t) and Ω 2 (t) :

Ω 1 (t) = x ∈ R 2 + / c(t -T ) < r * (x) ≤ ct Ω 2 (t) = x ∈ R 2 + / r * (x) ≤ c(t -T ) (13)
These two sets are represented on figures 6 and 7 for two values of t. Note that Ω 1 (t) is not empty as soon as t > (h + 2L)/c while Ω 2 (t) is not empty as soon as t > (h + 2L)/c + T .

(h + 2L, 0) ct c(t -T ) Ω 1 (t) Fig. 6. The sets Ω 1 (t) and Ω 2 (t), if (h + 2L)/c ≤ t < T + (h + 2L)/c. Ω 1 (t) Ω 2 (t) (h + 2L, 0) Fig. 7. The sets Ω 1 (t) and Ω 2 (t), if t ≥ T + (h + 2L)/c.
According to [START_REF] Lassas | On the existence and convergence of the solution of pml equations[END_REF], in order to derive an L ∞ estimate of e σ,L (., t), we need an upper bound for the quantity:

sup x ∈ Ω 1 (t) G N r (x, .) L 1 ( r * c , t) when t > h + 2L c ,
and for the quantity:

sup x ∈ Ω 2 (t) G N r (x, .) L 1 (t-T, t) when t > h + 2L c + T.
We thus first have to estimate L 1 -norms in time of G σ,L r (x, .). Estimating the cosine function by 1, we get :

|G σ,L r (x, t)| < e -A * (x,t) 2π t 2 -r * (x) 2 c 2 .
Let us remark that the function t → e -A * (x,t) is decreasing for t > 0, thus, recalling that, for x 1 < 0, Σ * (x 1 ) = 2Σ(L) = 2Lσ :

G σ,L r (x, .) L 1 ( r * c ,t) ≤ 1 2π e -A * (x, r * c ) t r * c dτ t 2 -r * 2 c 2 , = 1 2π e -2L c σ| cos θ * | Log ct + √ c 2 t 2 -r * 2 r * (14)
in the same way :

G σ,L r (x, .) L 1 (t-T,t) ≤ e -2Lσ| cos θ * | t-T r * 2π Log   ct + √ c 2 t 2 -r * 2 c(t -T ) + c 2 (t -T ) 2 -r * 2   . ( 15
) Since cos θ * (x) = (x 1 -2L -h)/r * (x) we easily verify that inf x ∈ Ω 1 (t) | cos θ * (x)| = 2L + h ct and inf x ∈ Ω 2 (t) | cos θ * (x)| r * (x) = 2L + h c 2 (t -T ) 2 , hence sup x ∈ Ω 1 (t) e -2L c σ| cos θ * (x) | = e -2L(2L+h)σ c 2 t
and sup

x ∈ Ω 2 (t) e -2Lσ sin θ * t-T r * = e -2L(2L+h)σ c 2 t
.

On the other hand, using the fact that the two functions:

r → ct r + c 2 t 2 r 2 -1 , r ∈ [0, ct], r → ct + √ c 2 t 2 -r 2 c(t -T ) + c 2 (t -T ) 2 -r 2 , r ∈ [0, c(t -T )], (t > T )
are respectively decreasing and increasing, we deduce that:

sup x ∈ Ω 1 (t) Log ct r * + c 2 t 2 r * 2 -1 = Φ(t)
and sup

x ∈ Ω 2 (t) Log ct + √ c 2 t 2 -r * 2 c(t -T ) + c 2 (t -T ) 2 -r * 2 = Φ(t).
It is now obvious to obtain :

             sup x ∈ Ω 1 (t) G σ,L r (x, .) L 1 ( r * c , t) ≤ 1 2π e -2L(2L+h)σ c 2 t Φ(t) sup x ∈ Ω 2 (t) G σ,L r (x, .) L 1 (t-T, t) ≤ 1 2π e - 2L(2L+h)σ c 2 (t-T ) Φ(t) (16)
Since, for t > 2L+h c + T , exp(-1/t) > exp(-1/(t -T )), it is easy to obtain the uniform estimates from the first inequality.

Numerical Results

The Green's function. The main difficulty, if one wants to represent numerically the function G σ,L r (x, .) is the singularity on the circle r * (x) = ct. To overcome this problem we introduce the relative error defined by :

γ σ,L r (x, t) = G σ,L r (x, .) G 0,L r (x, .) = e -A(x,t) cos [B(x, t)] . (17)
Note that G 0,L r (x, .) represents the field we would have obtained by imposing a Neumann condition on the line x 1 = L without adding a damping term. That is why γ σ,L is called a relative error. We choose h = 1 and c = 1 The figure 8 represents the level curves of γ σ,L for three different values of σ (σ=10, 20 and 50 from the left to the right). We remark that the amplitude of the error decrease strongly with σ : the error level is 0.6 for σ = 10, 0.3 for σ = 20 and 0.06 for σ = 50. When σ increases, the amplitude of the error concentrates more and more at the neighborhood of the layer. Moreover, its dependence with respect to the space variable is much more complicated because of the influence of the cosine function. The case of a source term. We have implemented a Matlab code to com-pute the convolution integral [START_REF] Collino | Optimizing the Perfectly Matched Layer[END_REF]. To validate our solution we have compared it to the one obtained by a finite difference code (based on a standard central finite difference scheme of leap-frog type, second order accurate in space and time).

In our experiment, the source function is a truncated first derivative of a Gaussian: On figures 11 and 12 we have compared both solutions at point (0.9, 0.1) as functions of time. The solid curves represent the "analytical" solution and the dashed curves the numerical one for two values of σ: σ = 1.5 and 3. As before the left pictures represent the total field while the right pictures represent the reflected field.

f (t) = d dt e -2πf 0 (t-t 0 ) 2 H(2t 0 -t), f 0 = 10, t 0 = 1/f 0 (18
L ∞ error estimates. On figure 13 we have compared the L ∞ norm of the reflected field (the solid curves) to the uniform estimates [START_REF] Collino | The Perfectly Matched Layer in curvilinear coordinates[END_REF] given by Theorem 4 (the dashed curves) for σ = 1, 10, 100. The source is a step function in time : The case of a rectangular domain surrounded by PML. We consider a domain [0 ; 1] 2 surrounded by a PML of width L = 0.1. We use a point source in space located at the center of the domain, the expression of the source in time is given by [START_REF] De Hoop | Analytic time-domain performance analysis of absorbing boundary conditions and perfectly matched layers[END_REF]. As for the case of the half-plane, we have compared our results to a numerical solution obtained by the same code as for the previous paragraph . Theses results are represented in figures 14, 15 and 16, respectively at time 7, 12 and 17 (take care to the fact that the color scale differs from one picture to another). Here again the results show a very good agreement between the analytical and numerical solution, even for long time, after several reflections. The Bérenger's PML were adapted to the aeroacoustics equations for instance by Hu [START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF] in 1996. Let us first recall his model, considering an horizontal uniform mean flow M. We have to rewrite the aeroacoustics equations whose unknowns are the pressure and the velocity (p, U = (u, v)),

f (t) = 1 if 0 ≤ t ≤ 2
                     ∂p ∂t + M ∂p ∂x 1 + ∂u ∂x 1 + ∂v ∂x 2 = 0, ∂u ∂t + M ∂u ∂x 1 + ∂p ∂x 1 = 0, ∂v ∂t + M ∂v ∂x 1 + ∂p ∂x 2 = 0, (19)
in a so-called "split form" with the unknowns (p 1 , p 2 , u, v 1 , v 2 ) where p 1 and p 2 (resp. v 1 and v 2 ) are non physical variables, whose sum is equal to the pressure (resp. the second component of the velocity). One then obtains the PML model by adding a zero-order absorption term, for instance in the equations containing the derivatives with respect to x 1 for a vertical layer :

                                           ∂p 1 ∂t + σ(x 1 )p 1 + M ∂p ∂x 1 + ∂u ∂x 1 = 0, ∂p 2 ∂t + ∂v ∂x 2 = 0, ∂u ∂t + σ(x 1 )u + M ∂u ∂x 1 + ∂p ∂x 1 = 0, ∂v 1 ∂t + σ(x 1 )v 1 + M ∂v ∂x 1 = 0, ∂v 2 ∂t + ∂p ∂x 2 = 0. ( 20
)
Though this problem is (weakly) well-posed (see [START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF]), it is unstable as soon as |M| > 0 : it may admit solution growing exponentially with the time variable.

Analysis of the instabilities via the slowness curve

The instabilities can be analyzed thanks to a plane wave analysis and to the slowness curves : we search for (ω,

k 1 , k 2 ) such that [p 1 , p 2 , u, v 1 , v 2 ] t (x 1 , x 2 , t) = [p 0 1 , p 0 2 , u 0 , v 0 1 , v 0 2 ] t e i(ωt+k 1 x 1 +k 2 x 2 ) .
It is well-known that, if σ = 0 (i.e for the aeroacoustics equations), (ω, k 1 , k 2 ) is solution of one of the two following dispersion equations :

• the dispersion equation of the advective waves :

(ω + Mk 1 ) 2 -k 2 1 -k 2 2 = 0; (21)
• the dispersion equation of the vorticity waves :

ω + Mk 1 = 0. ( 22
)
The slowness curves are the representation in the so-called slowness plane (k 1 /ω, k 2 /ω) of the solutions of these equations. One can easily show that the slowness curve associated to the advective waves is an ellipse whose center is (M/(1 -M 2 ), 0) and that the slowness curve associated to the vorticity is a line of equation k 1 /ω = -1/M. Let us also define the slowness vector and the group velocity respectively by : l = [k 1 /ω, k 2 /ω] t and v g = ∇ k ω. It is immediate to deduce the slowness vector from the slowness curve. Moreover the group velocity vector is orthogonal to the slowness curves (see figure 17). It is well-known that, as soon as |M| > 0, it exists waves whose first component of slowness and of group velocity vector have opposite signs (in red in figure 17).

It has been shown [START_REF] Hesthaven | On the analysis and construction of perfectly matched layers for the linearized Euler equations[END_REF][START_REF] Tam | Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains[END_REF][START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF] that the high-frequency instabilities are due to these so-called back-propagating modes. It is thus easy to deduce the set of unstable waves : the advective waves such that 0 < k 1 /ω < M/(1 -M 2 ) (see figure 18). An illustration of this instabilities is given in figure 19 for M = 0.5. As for the acoustics equation, it is possible to obtain an analytic expression of the fundamental solution of problem 20. Let us first remark that this problem can be formally rewritten as (the calculations are left to the reader) :

k 2 ω k 1 ω M 1-M 2 - 1 M v g v g v g v g O l l l l
                     ∂p ∂t + MD σ x 1 p + D σ x 1 u + ∂v ∂x 2 = 0, ∂u ∂t + MD σ x 1 u + D σ x 1 p = 0, ∂v ∂t + MD σ x 1 v + ∂p ∂x 2 = 0. (23)
Then, straightforward calculations lead us to the following decoupled system :

                     ∂ 2 p ∂t 2 + 2MD σ x 1 ∂p ∂t -(1 -M 2 )(D σ x 1 ) 2 p - ∂ 2 p ∂x 2 2 = 0, ∂u ∂t + MD σ x 1 u + D σ x 1 p = 0, ∂v ∂t + MD σ x 1 v + ∂p ∂x 2 = 0. ( 24
)
The first equation is nothing but the advective wave equation where the derivation with respect to x 1 has been replaced by the operator D σ x 1 . This equation, which we choose to study here, is the cause of the instabilities.

As for the acoustics case we consider the problem

             Find u : IR 2 × IR + → IR, null for t < 0 such that ∂ 2 p ∂t 2 +2M ∂ ∂t D σ x 1 p -(1 -M 2 )(D σ x 1 ) 2 p+ ∂ 2 ∂x 2 2 p = δ(x -x S )δ(t) in IR 2 × IR + . ( 25
) where x S = (-h, 0), σ(x 1 ) = 0 if x 1 ≤ 0 and σ(x 1 ) = 0 if x 1 > 0.
We define the functions r(x) and θ(x) ("generalized polar coordinates") as :

r 2 (x) = (x 1 + h) 2 (1 -M 2 ) 2 + x 2 2 1 -M 2 , cos θ(x) = x 1 + h r(1 -M 2 ) and sin θ(x) = x 2 r √ 1 -M 2 .
Note that the level lines of r(x) and θ(x) are ellipses and half-lines (see Fig. 20). We also introduce the two functions A(x, t) and B(x, t) defined by

r(x) = cst θ(x) = cst
A(x, t) = Σ(x 1 ) t r cos θ -M sin 2 θ and B(x, t) = Σ(x 1 ) sin θ t r + M cos θ 2 -1.
Theorem 6 The expression of the solution p(x, t) = G σ i (x, t) of problem (25) is given by :

G σ i (x, t) = H(t -r(1 -M cos θ)) 2π √ 1 -M 2 (t + Mr cos θ) 2 -r 2 , x 1 ≤ 0 G σ i (x, t) = H(t -r(1 -M cos θ)) 2π √ 1 -M 2 (t + Mr cos θ) 2 -r 2 e -A(x,t) cos [B(x, t)] , x 1 > 0. ( 26
)
For x 1 > 0 and t > t 0 = r(1 -M cos θ) we define the function

γ(x, t) = G σ i (x, t) G 0 i (x, t) = e -A(x,t) cos [B(x, t)] .
Theorem 7 For fixed x, the function t → γ(x, t) is bounded : γ(x, t) ≤ C(x). However, given a direction θ 0 , let us denote

r max (t) = t (1 -M cos θ 0 ) , then, γ(r max (t), θ 0 , t) > e α(θ 0 )t with α(θ 0 ) > 0 if cos θ 0 < M and α(θ 0 ) ≤ 0 if cos θ 0 ≥ M. Consequently the L ∞ (IR + )-norm of γ(x, t) grows exponentially with t : γ(., t) L ∞ ≥ Ce αt , α > 0.
The kind of instability expressed by this theorem is often called "convective instability" because the instability is "moving with the wave front". To illustrate this point we have represented the variations of γ with respect to r for a given unstable direction θ 0 and for four values of t in figure 21.

We have also implemented a Matlab code to compute the solution of problem [START_REF] Hesthaven | On the analysis and construction of perfectly matched layers for the linearized Euler equations[END_REF] where the source function is once again a truncated first derivated of a Gaussian. In the following experiment the frequency of the source is 1 and the source is located at point (0, -2). Figures [START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF] and ( 23) represents the propagation of an advective wave, respectively with M = 0.5 and M = 0.9. We have represented in red the half lines θ(x) = arccos(M) and θ(x) = -arccos(M) that define the set of unstable points.

Proof of the theorems

Expression of the analytic solution. Let us first apply to G σ i the Laplace transform in t and the Fourier transform in x 2 . The resulting function x 1 → G σ i (x 1 , k, s) satisfies the ordinary differential equation :

-(1 -M 2 ) s s + σ d dx 1 s s + σ d G σ i dx 1 + 2Ms s s + σ ∂ G σ i ∂x 1 + (k 2 + s 2 ) G σ i = δ(x 1 + h) (27) We denote λ(x 1 ) =          M + 1 + (1 -M 2 ) k 2 s 2 , if x 1 < -h, M -1 + (1 -M 2 ) k 2 s 2 , if x 1 > -h, G σ i (x 1 , k, s) has the following form G σ i (x 1 , k, s) = A(k, s) e s(x 1 +h)+Σ(x 1 ) 1-M 2 λ(x 1 )
Taking into account the jump condition

-(1 -M 2 ) s s + σ 2 ∂ G σ i ∂x 1 + 2M s s + σ G σ i x 1 =-h = 1, ( 28
) one obtains, if σ(-h) = 0, A(k, s) = 1 2 s 2 + (1 -M 2 )k 2 . If x 1 is negative, Σ(x 1 ) = 0 and G σ i (x 1 , k, s) = e s(x 1 +h) 1-M 2 λ(x 1 ) 2 1 + (1 -M 2 ) k 2 s 2
. This is nothing but the Fourier-Laplace transform of the fundamental solution of the advective wave equation :

G σ i (x, t) = H(c t -r(1 -M cos θ)) 2π (t + Mr cos θ) 2 -r 2 , x 1 ≤ 0.
Let us now consider

x 1 > 0 so that λ(x 1 ) = M -1 + (1 -M 2 ) k 2 s 2 and
G σ i (x 1 , k, s) = e s(x 1 +h)+Σ(x 1 ) 1-M 2 M -1+(1-M 2 ) k 2 s 2 2 s 2 + (1 -M 2 )k 2
Applying the x 2 inverse Fourier transform to G σ i (x 1 , k, s), we get

G σ i (x 1 , x 2 , s) = 1 4π +∞ -∞ e s(x 1 +h)+Σ(x 1 ) 1-M 2 M -1+(1-M 2 ) k 2 s 2 s 2 + (1 -M 2 )k 2 e -ikx 2 dk, (29)
We apply the change of variable k = ps √ 1-M 2 and obtain:

G σ i (x 1 , x 2 , s) = 1 4π +∞ -∞ e Σ(x 1 ) M - √ 1+p 2 1-M 2 e s ( M - √ 1+p 2 ) 1-M 2 (x 1 +h)-ip x 2 √ 1-M 2 (1 -M 2 )(1 + p 2 ) dp. ( 30
)
We now have to search a path in the complex plane such that

M - √ 1 + p 2 1 -M 2 (x 1 + h) -ip x 2 √ 1 -M 2 = -t, for t > 0, (31)
For t ≥ t 0 = r (1 -M cos θ) we define the functions γ ± (t) by

γ ± (t) ≡ -i t r + M cos θ sin θ ± cos θ ( t r + M cos θ) 2 -1,
and we introduce the Cagniard-de Hoop contour Γ defined as:

Γ = Γ + ∪ Γ -, Γ ± = p = γ ± (t) . ( 32
)
The contour Γ is very similar to the one exposed at the previous section (consider denoting T = t + Mr cos θ) : in particular it does not intercept the branch cut of the function (1+p 2 ) -1 2 . Then, as for the acoustics case, we easily deduce :

G σ i (x 1 , x 2 , s) = - 1 4π Γ e Σ(x 1 ) M - √ 1+p 2 1-M 2 e s ( M - √ 1+p 2 ) 1-M 2 (x 1 +h)-ip x 2 √ 1-M 2 √ 1 + p 2 dp.
We use the parameterizations p = γ + (t) and p = γ -(t), for t ≥ t 0 , respectively along Γ + and Γ -and remark that:

• M - √ 1 + p 2 1 -M 2 (x 1 + h) -ip x 2 √
1 -M 2 = -t, (by construction),

• dp

(1 + p 2 ) 1 2 = ± dt (t + Mr cos θ) 2 -r 2 1 2 , on Γ ± , • Σ(x 1 ) M - √ 1 + p 2 1 -M 2 = -A(x, t) ∓ iB(x, t), on Γ ±
Therefore, since t goes from +∞ to t 0 on Γ + and from t 0 to +∞ on Γ -:

G σ i (x 1 , x 2 , s) = 1 4π +∞ t 0 e -A(x,t) cos (B(x, t)) √ 1 -M 2 (t + Mr sin θ) 2 -r 2 e -st dt.
It is then easy to conclude.

Properties of the function γ. For fixed x, the function B(x, t) is bounded by 1 in time. Moreover the function t → A(x, t) is obviously increasing with t on [t 0 ; +∞] and its minimum is then A(x, t 0 ) = Σ(x 1 )(cos θ -M). Then max t≥t 0 γ(x, t) ≤ e -Σ(x 1 )(cos θ-M ) .

For given t we denote by Ω(t) the support of γ(., t), it is a disk of center (Mt -h, 0) and of radius t. Let us remark that the point of generalized polar coordinates (r max (t), θ 0 ) belongs to ∂Ω t . Denoting

x max 1 (t) = (1 -M 2 )r max (t) cos θ 0 -h = max ∆(θ 0 ,t) x 1 ,
it is easy to check that

γ(r max (t), θ 0 , t) = e -Σ(x max 1 (t))(cos θ 0 -M ) ,
since B(r max (t), θ 0 , t) = 0, then, since x max 1 (t) grows linearly with t and Σ grows with x 1 it is easy to conclude the proof 3 Stabilized PML for advective acoustics Various solutions have been proposed in the literature to prevent instabilities in the PML. Hu [START_REF] Hu | On absorbing boundary conditions for linearized euler equations by a perfectly matched layer[END_REF] used a low-pass filter inside the absorbing layer. Tam, Auriault and Cambuli [START_REF] Tam | Perfectly matched layer as an absorbing boundary condition for the linearized Euler equations in open and ducted domains[END_REF] proposed to use selective damping coefficients. Hesthaven [START_REF] Hesthaven | On the analysis and construction of perfectly matched layers for the linearized Euler equations[END_REF] decelerated progressively the flow in the layer. Lions, Metral and Vacus [START_REF] Lions | Well-posed absorbing layer for hyperbolic problems[END_REF] defined a new layer by regularizing the damping terms. However all these methods may be not perfectly matched and, although these problems are well-posed, their stability has not be proven. Recently a significative advance has been achieved by Abarbanel, Gottlieb and Hesthaven [START_REF] Abarbanel | Well-posed perfectly matched layers for advective acoustics[END_REF] and by Hu [START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF] who has presented a stable and perfectly matched model, close to the one we are going to present (though the two models have been obtained independently). Hagström [START_REF] Hagström | Absorbing layers and radiation boundary condition for jet flow simulations[END_REF], using a different method, has also constructed a similar layer.

The advective wave equation

We are first concerned by stabilizing the PML for the advective equation whose associated slowness curve is the ellipse represented in the figure 18. As we said before, the instabilities produced by the PML can be explained by the fact that the center of this ellipse is not the origin of the axes. Therefore an idea for stabilizing the PML consists in using a change of variable to translate the ellipse by the vector (-M/(1 -M 2 ), 0) for its center to be located at the origin (see figure 24). Therefore we introduce the new unknowns k ⋆ 1 and ω

⋆ k 2 ω k 1 ω M 1 -M 2
such that k ⋆ 1 ω ⋆ = k 1 ω - M 1 -M 2 . If we set ω ⋆ = ω this relation becomes k ⋆ 1 = k 1 - M 1 -M 2 ω.
This change of variable can also be expressed through the time-space variables :

∂ ∂x 1 - M 1 -M 2 ∂ ∂t =⇒ ∂ ∂x ⋆ 1 ,
which is obtained by using the new variables :

x ⋆ 1 = x 1 , x ⋆ 2 = x 2 , and t ⋆ = t - M 1 -M 2 x 1 . (33)
Note that this change of variable does not affect the spatial coordinates and consequently the interface between the computational domain and the layer. The choice ω ⋆ = ω has been made to obtain this property.

Construction of the stabilized PML model. Let us first show how to obtain an infinite PML in IR 2 + (x 1 > 0). The steps of the methods are :

(1) Rewrite the problem as a transmission problem between IR 2 -and IR 2 + :

∂ 2 p ∂t 2 + 2M ∂p ∂t∂x 1 -(1 -M 2 ) ∂ 2 p ∂x 2 1 - ∂ 2 p ∂x 2 2 = δ(x -x s )δ(t), x 1 < 0, (34a)
∂ 2 p ∂t 2 + 2M ∂p ∂t∂x 1 -(1 -M 2 ) ∂ 2 p ∂x 2 1 - ∂ 2 p ∂x 2 2 = 0, x 1 > 0, (34b) p(0 -, x 2 , t) = p(0 + , x 2 , t), (34c)
∂p ∂x 1 (0 -, x 2 , t) = ∂p ∂x 1 (0 + , x 2 , t). ( 34d
)                                   
(2) Apply the change of variable

p ⋆ (x 1 , x 2 , t) = p(x 1 , x 2 , t - M 1 -M 2 x 1 ),
in the PML to obtain the new transmission problem :

∂ 2 p ∂t 2 + 2M ∂p ∂t∂x 1 -(1 -M 2 ) ∂ 2 p ∂x 2 1 - ∂ 2 p ∂x 2 2 = δ(x -x s )δ(t), x 1 < 0, (35a) 1 1 -M 2 ∂ 2 p ⋆ ∂t 2 -(1 -M 2 ) ∂ 2 p ⋆ ∂x 2 1 - ∂ 2 p ⋆ ∂x 2 2 = 0 x 1 > 0, (35b) p(0, x 2 , t) = p ⋆ (0, x 2 , t), (35c)
∂p ∂x 1 (0, x 2 , t) = ∂p ⋆ ∂x 1 (0, x 2 , t) + M 1 -M 2 ∂p ⋆ ∂t (0, x 2 , t). ( 35d
)                                   
Note that the advective wave equation has been turned into an equation very close to the classical wave equation by this change of variable.

Let us now define the function ρ(x) by ρ(x) = 1 if x 1 < 0 and ρ(x) = 1/(1 -M 2 ) otherwise.

Theorem 8

The following energy is conserved :

E(t) = 1 2 IR 2   ρ(x) ∂p ∂t 2 +(1 -M 2 ) ∂p ∂x 1 2 + ∂p ∂x 2 2   .
Proof. Let us respectively multiply (35a) and (35b) by ∂p/∂t and ∂p ⋆ /∂t. After having integrated both this equations over their respective defini-tion domain (IR 2 -and IR 2 + ) and added them, we check that :

dE dt (t) + IR 2 - 2M ∂p ∂t∂x 1 ∂p ∂t + Γ (1 -M 2 ) ∂p ⋆ ∂x 1 ∂p ⋆ ∂t - ∂p ∂x 1 ∂p ∂t + Γ ∂p ⋆ ∂x 2 ∂p ⋆ ∂t - ∂p ∂x 2 ∂p ∂t = 0
where Γ is the line of equation x 1 = 0. Using the continuity of p (and consequently of its derivatives with respect to x 2 ant t) it is obvious that the last integral in this equation is zero. Furthermore

IR 2 - M ∂p ∂t∂x 1 ∂p ∂t = Γ M ∂p ∂t ∂p ∂t - IR 2 - M ∂p ∂t ∂p ∂t∂x 1 .
Then

IR 2 - 2M ∂p ∂t∂x 1 ∂p ∂t + Γ (1 -M 2 ) ∂p ⋆ ∂x 1 ∂p ⋆ ∂t - ∂p ∂x 1 ∂p ∂t = Γ ∂p ∂t M ∂p ∂t + (1 -M 2 ) ∂p ⋆ ∂x 1 - ∂p ∂x 1
We deduce from (35d) that this expression is zero :

dE dt (t) = 0.
(3) It only remains to replace formally the derivatives with respect to x 1 by the operator D σ x 1 :

∂ 2 p ∂t 2 + 2M ∂p ∂t∂x 1 -(1 -M 2 ) ∂ 2 p ∂x 2 1 - ∂ 2 p ∂x 2 2 = δ(x -x s ), x 1 < 0, (36a) 1 1 -M 2 ∂ 2 p ⋆ ∂t 2 -(1 -M 2 )(D σ x 1 ) 2 p ⋆ - ∂ 2 p ⋆ ∂x 2 2 = 0 x 1 > 0, (36b) p(0, x 2 , t) = p ⋆ (0, x 2 , t), (36c)
∂p ∂x 1 = D σ x 1 p ⋆ (0, x 2 , t) + M 1 -M 2 ∂p ⋆ ∂t . ( 36d
)                                   
Analysis of the new model. Let us define the functions r(x) and θ(x) by

r 2 (x) = (x 1 + h) 2 (1 -M 2 ) 2 + x 2 2 1 -M 2 , cos θ(x) = x 1 + h r(1 -M 2 ) and sin θ(x) = x 2 r √ 1 -M 2 .
Let A(x, t) and B(x, t) be defined by

A(x, t) = Σ(x 1 ) t + Mh 1 -M 2 cos θ r > 0 and B(x, t) = Σ(x 1 ) sin θ r t + Mh 1 -M 2 2 -1.
Theorem 9 The expression of the solution (p, p ⋆ )(x, t) = (G i , G ⋆ i )(x, t) of problem (36) is given by :

G i (x, t) = H(t -r(1 -M cos θ)) 2π √ 1 -M 2 (t + Mr cos θ) 2 -r 2 , x 1 ≤ 0 G ⋆ i (x, t) = H(t + M h 1-M 2 -r) 2π √ 1 -M 2 t + M h 1-M 2 2 -r 2 e -A(x,t) cos [B(x, t)] , x 1 > 0. ( 37
)
This expression is similar to one exposed at the previous section. The différence lies in the fact that the function A is always positive which guarantee the stability of the layer.

Proof via the Cagniard-de Hoop technique. After applying a Laplace transform in t and a Fourier transform in x 2 we obtain the following differential system :

(s 2 + k 2 ) G i + 2sM d G i dx 1 -(1 -M 2 ) d 2 G i dx 1 = δ(x 1 + h), x 1 < 0,(38a) s 2 1 -M 2 + k 2 G ⋆ i -(1 -M 2 ) s s + σ d dx 1 s s + σ d G i dx 1 = 0 x 1 > 0,(38b) G i (0) = G ⋆ i (0), (38c)
d G i dx 1 (0) = s s + σ d G ⋆ i dx 1 (0) + Ms 1 -M 2 G ⋆ i (0). ( 38d
)                                     
After calculations (similar to the ones of section 2) we obtain :

G i (x 1 , k, s) = e s(x 1 +h) 1-M 2 λ(x 1 ) 2 s 2 + (1 -M 2 )k 2 , G ⋆ i (x 1 , k, s) = e 1 1-M 2 shM -x 1 +h+ Σ(x 1 ) s √ s 2 +(1-M 2 )k 2 2 s 2 + (1 -M 2 )k 2 ,
where λ has been defined at the section 2. As G i is once again the Fourier-Laplace transform of the classical advective wave equation it is obvious that the new PML does not produce parasitical reflection. Applying the x 2 inverse Fourier transform to G ⋆ i (x 1 , k, s), we get

G ⋆ i (x 1 , x 2 , s) = 1 4π +∞ -∞ e 1 1-M 2 shM -x 1 +h+ Σ(x 1 ) s √ s 2 +(1-M 2 )k 2 2 s 2 + (1 -M 2 )k 2 e -ikx 2 dk, (39)
G ⋆ i (x 1 , x 2 , s) = 1 4π +∞ -∞ e -Σ(x 1 ) √ 1+p 2 1-M 2 e s ( M h-(x 1 +h) √ 1+p 2 ) 1-M 2 -ip x 2 √ 1-M 2 (1 -M 2 )(1 + p 2 ) dp. ( 40
)
We now have to search a path in the complex plane such that

Mh -(x 1 + h) √ 1 + p 2 1 -M 2 -ip x 2 √ 1 -M 2 = -t, for t > 0, (41)
For t ≥ t 0 = r -M h 1-M 2 > 0 we define the functions γ ± (t) by

γ ± (t) ≡ -i t r + Mh r(1 -M 2 ) sin θ ± cos θ t r + Mh r(1 -M 2 ) 2 -1,
and we introduce the Cagniard-de Hoop contour Γ defined as:

Γ = Γ + ∪ Γ -, Γ ± = p = γ ± (t) . (42)
As for the acoustics case we easily deduce :

G ⋆ i (x 1 , x 2 , s) = - 1 4π Γ e -Σ(x 1 ) √ 1+p 2 1-M 2 e s ( M h-(x 1 +h) √ 1+p 2 ) 1-M 2 -ip x 2 √ 1-M 2 (1 -M 2 )(1 + p 2 ) dp.
We use the parameterizations p = γ + (t) and p = γ -(t), for t ≥ t 0 , respectively along Γ + and Γ -and remark that:

• Mh -(x 1 + h) √ 1 + p 2 1 -M 2 -ip x 2 √ 1 -M 2 = -t, ( by
construction),

• dp

(1 + p 2 ) 1 2 = ± dt t + M h 1-M 2 2 -r 2 1 2 , on Γ ± , • -Σ(x 1 ) √ 1 + p 2 1 -M 2 = -A(x, t) ∓ iB(x, t), on Γ ±
Therefore, since t goes from +∞ to t 0 on Γ + and from t 0 to +∞ on Γ -:

G i (x 1 , x 2 , s) = 1 4π +∞ t 0 e -A(x,t) cos (B(x, t)) √ 1 -M 2 t + M h 1-M 2 2 -r 2 e -st dt.
It is then easy to conclude.

Numerical illustration. In figure 25 we compare three numerical experiments :

(1) (top picture) the propagation of an advective wave in a large computational domain : (2) (middle picture) the same experiment but the change of variable (33) has been applied in the left and right hand layers ; (3) (bottom picture) in the last one a damping term has been added in the layers. The last experiment shows that the waves are properly damped in the PML.

Extension to the linearized Euler equations

Looking at the shape of slowness curve associated to the vorticity waves (a line orthogonal to the axis k 1 /ω), we can reasonably assume that the change of variable (33), which translate the slowness curves of M/(1 -M 2 ) to the right (see figure 26), does not generate instabilities.

k y ω k x ω - 1 M M 1 -M 2
∂s ∂t + (1 + M) ∂s ∂x 1 + ∂v ∂x 2 = 0, (43a)
∂d ∂t -(1 -M) ∂d ∂x 1 + ∂v ∂x 2 = 0, (43b)
∂v ∂t + M ∂v ∂x 1 + 1 2 ∂s ∂x 2 + ∂d ∂x 2 = 0. ( 43c
)                           
Then the three steps of the construction of a stable PML are :

(1) Write a transmission problem between x 1 < 0 and x 1 > 0.

(2) Apply the change of variable (33) to obtain the following new equations in the domain x 1 > 0 :

1 1 -M ∂s ⋆ ∂t + (1 + M) ∂s ⋆ ∂x 1 + ∂v ⋆ ∂x 2 = 0, (44a) 1 1 + M ∂d ⋆ ∂t -(1 -M) ∂d ⋆ ∂x 1 + ∂v ⋆ ∂x 2 = 0, (44b) 1 1 -M 2 ∂v ⋆ ∂t + M ∂v ⋆ ∂x 1 + 1 2 ∂s ⋆ ∂x 2 + ∂d ⋆ ∂x 2 = 0, ( 44c
)                         
that are coupled to (43) by the continuity conditions :

s ⋆ (0, x 2 , t) = s(0, x 2 , t), d ⋆ (0, x 2 , t) = d(0, x 2 , t), v ⋆ (0, x 2 , t) = v(0, x 2 , t).
(3) Apply to (44) the usual PML model.

Remark 10

The well-posedness of the transmission problem obtained after step 2 is guaranteed by the following theorem, whose proof, omitted here, is similar to the one of theorem [START_REF] Collino | The Perfectly Matched Layer in curvilinear coordinates[END_REF].

Theorem 11

The following energy is conserved :

E 2 (t) = 1 2 IR- s 2 + d 2 + 2v 2 + 1 2 IR + s ⋆2 1 -M + d ⋆2 1 + M + 2 v ⋆2 1 -M 2 . (45)
Comparison with the Hu's model. The PML recently proposed by Hu are close to the one we have presented here. Both are based on the same change of variable, the difference lies in the fact that Hu applies this change of variable in the whole domain whereas we only apply it for x 1 > 0. Therefore he does not need to impose transmission condition at the interface. Then he replace formally the derivatives with respect to x 1 by the operator D σ x 1 . Finally he use the inverse change of variable in the whole domain to come back to the original variables. This inverse change of variable makes the operator D σ x 1 slightly more complicated.

Numerical Results

In figure (

A.1 Construction of a vertical PML for the advective equation

A non-parallel mean flow is characterized by the Mach vector M = [M 1 ; M 2 ] t and we set M = M . The advective wave equation is then : (1) Write a transmission problem between x 1 < 0 and x 1 > 0.

∂ 2 p ∂t 2 +2M 1 ∂ 2 p ∂t∂x 1 +2M 2 ∂ 2 p ∂t∂x 2 +2M 1 M 2 ∂ 2 p ∂x 1 ∂x 2 -(1-M 2 1 ) ∂ 2 p ∂x 2 1 -(1-M 2
(2.1) In the region x 1 > 0, get rid of the convective term by applying a change of variable :

p ⋆ (x 1 , x 2 , t) = p x 1 , x 2 , t - M 1 1 -M 2 x 1 - M 2 1 -M 2 x 2 .
The new slowness curve (see figure A.2) still gives rise to instability if one applies classical PML. This is due to the presence of the cross-derivative with respect to x 1 and x 2 .

(2.2) In the region x 1 > 0 we apply the second change of variable (which let invariant the axis x 1 = 0 ) x1 = αx 1 and x2 = x 2 + βx 1 .

Writing the equations in the new variables, one observes that the cross derivatives disappear if α and β are related by :

β = - M 1 M 2 1 -M 2

A.2 The corner problem

The solution we have adopted for corners is much less justified from a numerical point of view but appears to give satisfactory results (see Fig. A.4). Let us for instance suppose that the computational domain is the quarter-plane IR -× IR -. In the horizontal and vertical layers, according to the previous section, we use the following models

1 1 -M 2 ∂ 2 p ∂t 2 -(1 -M 2 1 )(D σ x 1 ) 2 p -(1 - M 2 2 1 -M 2 1 ) ∂ 2 p ∂x 2 2 = 0 (vertical) 1 1 -M 2 ∂ 2 p ∂t 2 -(1 - M 2 1 1 -M 2 2 ) ∂ 2 p ∂x 2 1 -(1 -M 2
2 )(D σ x 2 ) 2 p = 0 (horizontal).

In the corner we have chosen to use the following equation in which variables x 1 and x 2 are not treated in a symmetric way! )(D σ x 2 ) 2 p = 0

We refer to [START_REF] Diaz | Approches analytiques et numériques de problèmes de transmission en propagation d'ondes en régime transitoire. application au couplage fluidestructure et aux méthodes de couches parfaitement adaptées[END_REF] and present in figure A.4 a numerical result obtained with this method.

Fig. 1 .

 1 Fig. 1. Computational domain and Perfectly Matched Layer

Fig. 8 .

 8 Fig. 8. x → γ σ,0.1 r

  )On figures 9 and 10 we have compared the "analytical" solution (top picture in each figure) to the numerical one (bottom picture in each figure) for two values of σ : σ = 1.5 and 3. In each picture we represent the level lines of the solution at time t = 0.4. The left pictures represent the total field while the right pictures represent the reflected field (the error). For the representation, the reflected field has been amplified by a factor which depends on σ : 4 for σ = 1.5 and 10 for σ = 3. In each case, the results reveal a very good agreement between the two solutions.

Fig. 9 .

 9 Fig. 9. Comparison between analytical (top picture) and numerical (bottom picture) solution σ = 1.5.

Fig. 10 .

 10 Fig. 10. Comparison between analytical (top picture) and numerical (bottom picture) solution σ = 3.

Fig. 11 .

 11 Fig. 11.a. The total field.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. Comparison between analytical (top picture) and numerical (bottom picture) solution, σ = 1.5

Fig. 12 .Fig. 13 Fig. 13 Fig. 13

 12131313 Fig. 12. Comparison between analytical (top picture) and numerical (bottom picture) solution, σ = 3 very sharp for σ = 1 and becomes less accurate (although quite acceptable) when σ increases.

Fig. 13 .

 13 Fig. 13. Error estimates (dashed curve) and effective error (solid curve)

Fig. 14 . 7 Fig. 15 . 12 Fig. 16 . 17 2

 14715121617 Fig. 14. The analytical (left picture) and numerical (right picture) solutions at time 7

Fig. 17

 17 Fig. 17. The slowness and group velocity vectors

Fig. 18 .Fig. 19 .

 1819 Fig. 17. The slowness and group velocity vectors

Fig. 20 .Fig. 21 .

 2021 Fig. 20. The level lines of r(x) and θ(x)

Fig. 22 . 5 Fig. 23 .

 22523 Fig. 22. Instabilities in the PML with M = 0.5

Fig. 24 .

 24 Fig.[START_REF] Bécache | On the analysis of Bérenger's perfectly matched layers for Maxwell's equations[END_REF]. The slowness curve before the change of variable (in blue/solid) and after (in magenta/dashed)

  We apply the change of variable k = ps √ 1-M 2 and obtain:

Fig. 25 .

 25 Fig. 25. Comparison between three experiments If we compare the two first experiments we remark that the change of variable (33) slows down the wave in the right hand layer whereas it accelerate it in the left hand one (actually, this change of variable "makes the flow vanish"). The last experiment shows that the waves are properly damped in the PML.

Fig. 26 .

 26 Fig. 26. The slowness curve before (in blue/solid) and after (in magenta/dashed) the change of variable It is useful to introduce the new unknowns s = p + u and d = p -u, and to rewrite the aeroacoustics equation as :

  27) we represent the propagation of an advective wave in an uniform mean flow with M = 0.5 in open domain. The domain of interest is [-5 ; 5] 2 and we impose an initial condition p = exp(-4(x 2 1 +x 2 2) ln(2)). The width of the PML is L = 1 and the profile of the damping coefficient is σ(x) = σ 0 x 2 /D. We have represented the level lines of the pressure at time t = 1, 3, 8, 30 and 50 We have repeated the same experiment but with an initial condition

Fig. 27 .

 27 Fig. 27. Propagation of an advective wave t = 1, 3, 8, 30 and 50u = x 1 exp(-(x 2 1 + x 2 2 ) ln(2)) and v = -x 2 exp(-(x 2 1 + x22 ) ln(2)) designed to simulate the propagation of a vorticity wave. Figure28represent the level lines of the modulus of the velocity at time t = 1, 3, 7, 11 and 50.

Fig. 28 .

 28 Fig. 28. Propagation of a vorticity wave t= 1, 3, 7, 11 and 50 A The case of a non-parallel mean-flow

2 )

 2 associated to this equation is represented in figure A.1 (the set of unstable points for a vertical layer has been highlighted in red), it is an ellipse of center C = (M 1 /(1 -M 2 ), M 2 /(1 -M 2 )). The construction of stabilized PML's is inspired by what has been done in the case of an horizontal mean flow. The novelty lies in the apparition of a new step (step (2.2))

1 α.M 2 1 -M 2 1 x 1 ( 3 )Fig. A. 3 .

 121133 Fig. A.1. The slowness curve associated to the advective waves

Fig. A. 4 .

 4 Fig. A.4. The case of a non-parallel mean flow, t= 1 3 7 15 and 50


Cécile Picard
email: cecile.picard@sophia.inria.fr

François Faure
email: francois.faure@inrialpes.fr

George Drettakis
email: george.drettakis@sophia.inria.fr

Paul G Kry
email: kry@cs.mcgill.ca

A ROBUST AND MULTI-SCALE MODAL ANALYSIS FOR SOUND SYNTHESIS

This paper presents a new approach to modal synthesis for rendering sounds of virtual objects. We propose a generic method for modal analysis that preserves sound variety across the surface of an object, at different scales of resolution and for a variety of complex geometries. The technique performs automatic voxelization of a surface model and automatic tuning of the parameters of hexahedral finite elements, based on the distribution of material in each cell. The voxelization is performed using a sparse regular grid embedding of the object, which easily permits the construction of plausible lower resolution approximations of the modal model. With our approach, we can compute the audible impulse response of a variety of objects. Our solution is robust and can handle non-manifold geometries that include both volumetric and surface parts, such as those used in games, training simulations, and other interactive virtual environment.

INTRODUCTION

Our goal is to realistically model sounding objects for animated real-time virtual environments. To achieve this, we propose a robust and flexible modal analysis approach that efficiently extracts modal parameters for plausible sound synthesis while also focusing on efficient memory usage.

Modal synthesis models the sound of an object as a combination of sinusoids, each of which oscillates independently of the others. Modal synthesis approaches are only accurate for sounds produced by linear phenomena, but they can compute these sounds in real-time. Modal synthesis requires the computation of a partial eigenvalue decomposition of the system matrices, which is relatively expensive. For this reason, modal analysis is performed in a preprocessing step. The eigenvalues and eigenvectors strongly depend on the geometry, material and scale of the sounding object. Therefore, modeling numerous sounding objects can rapidly become prohibitively expensive. In addition, this processing step can be subject to computation problems; in particular, when the geometries are non-manifold.

We propose a new approach to efficiently extract modal parameters for any given geometry, overcoming many of the afore mentioned limitations. Our method employs bounding voxels of a given shape at arbitrary resolution for hexahedral finite elements. The advantages of this technique are the automatic voxelization of a surface model and the automatic tuning of the finite element method (FEM) parameters based on the distribution of material in each cell. A particular advantage of this approach is that we can easily deal with non-manifold geometry which includes both volumetric and surface parts. These kinds of geometries cannot be processed with traditional approaches which use a tetrahedralization of the model (e.g., [START_REF] James | Synthesizing sounds from rigid-body simulations[END_REF]). Likewise, even with solid watertight geometries, complex details often lead to poorly shaped tetrahedra and numerical instabilities; in contrast, our approach does not suffer from this problem. Our specific contribution is the adaptation of the multi-resolution hexahedral embedding technique to modal analysis for sound synthesis. Most importantly, our solution preserves variety in what we call the Sound Map, that is, the changes in sound across the surface of the sounding object.

BACKGROUND

Related Work

The traditional approach to creating soundtracks for interactive physically based animations is to directly play-back pre-recorded samples, for instance, synchronized with the contacts reported from a rigid-body simulation. Due to memory constraints, the number of samples is limited, leading to repetitive audio. Moreover, matching sampled sounds to interactive animation is difficult and often leads to discrepancies between the simulated visuals and their accompanying soundtrack. Finally, this method requires each specific contact interaction to be associated with a corresponding pre-recorded sound, resulting in a time-consuming authoring process.

Work by Adrien [START_REF] Adrien | The missing link: modal synthesis[END_REF] describes how effective digital sound synthesis can be used to reconstruct the richness of natural sounds. There has been much work in Computer music [START_REF] Perry | Real Sound Synthesis for Interactive Applications[END_REF][START_REF] Iovino | Recent work around modalys and modal synthesis[END_REF] and computer graphics [START_REF] James | Synthesizing sounds from rigid-body simulations[END_REF][START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF][START_REF] Kees Van Den | Foley automatic: physically-based sound effects for interactive simulation and animation[END_REF] exploring methods for generating sound based on physical simulation. Most approaches target sounds emitted by vibrating solids. Physically based sounds require significantly more computation power than recorded sounds. Thus, brute-force sound simulation cannot be used for real-time sound synthesis. For interactive simulations, a widely used solution is to apply vibrational parameters obtained through modal analysis. Modal data can be obtained from simulations [START_REF] James | Synthesizing sounds from rigid-body simulations[END_REF][START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF] or extracted from recorded sounds of real objects [START_REF] Kees Van Den | Foley automatic: physically-based sound effects for interactive simulation and animation[END_REF]. The technique presented in this paper is more closely related to the work of O'Brien et al. [START_REF] James | Synthesizing sounds from rigid-body simulations[END_REF], which extends modal analysis to objects that are neither simple shapes nor available to be measured.

The computation time required by current methods to preprocess the modal analysis prevents it from being used for real-time rendering. Work of Maxwell and Bindel [START_REF] Maxwell | Modal parameter tracking for shape-changing geometric objects[END_REF] address interactive sound synthesis and how the change of the shape of a finite element model affects the sound emission. Concerning the performance of mode-based computations, a fast sound synthesis approach that exploits the inherent sparsity of modal sounds in the frequency domain has recently been introduced by Bonneel et al. [START_REF] Bonneel | Fast modal sounds with scalable frequency-domain synthesis[END_REF]. Our technique tackles computational efficiency by proposing a multi-scale resolution approach of the modal analysis, managing the amount of modal data according to the memory requirements.

Modal Analysis

Modal sound synthesis is a physically based synthesis approach that consists of solving the governing equations of motion of a sounding system. The solution of these equations becomes complicated when the size of the system is large or when the forcing functions of the system are non-periodic. A system of n degreesof-freedom is governed by a set of n coupled ordinary differential equations of second order. By expressing the deformation of the object as linear combinations of normal modes, the equations of motion are uncoupled and the solution for object vibration can be easily computed. In order to decouple the damped system into single degree-of-freedom oscillators, we assume Rayleigh damping (see, for instance, [START_REF] Bathe | Finite element procedures in engineering analysis[END_REF]). Using the finite element method, we obtain the general form of the system for eigendecomposition, from which the modal parameters, i.e., frequencies, dampings, and corresponding gains are extracted. For our approach, the calculations for modal parameters are similar to the ones presented in the paper of O'Brien et al. [START_REF] James | Synthesizing sounds from rigid-body simulations[END_REF] and we refer the reader to this work for additional information.

METHOD

In the case of small elastic deformations, rigid motion of an object does not interact with the objects's vibrations. On the other hand, we assume that small-amplitude elastic deformations will not significantly affect the rigid-body collisions between objects. For these reasons, the rigid-body behavior of the objects can be modelled in the same way as animation without audio generation.

Deformation Model

Our implementation uses the Sofa Framework1 for rigid-body simulation. This choice was motivated by the ease with which it could be extended for our purpose. The main feature of SOFA compared with other libraries is its high flexibility while maintaining efficiency. SOFA is an open-source C++ library for physical simulation. It can be used as an external library in another program, or using one of the associated GUI applications. It allows the use of multiple interacting geometrical models of the same object, typically, a mechanical model with mass and constitutive laws and a collision model with simple geometry. A visual model with detailed geometry and rendering parameters is also integrated, where each model can be designed independently of the others. During run-time, consistency is maintained using mappings. Additionally, ALGORITHM 1. Algorithm for modal parameters extraction.

1. Compute mass and stiffness at desired mechanical level 2. Assemble the mass and the stiffness matrices 3. Modal analysis: solve the eigenproblem 4. Store eigenvalues and eigenvectors for sound synthesis

The model approximates the motion of the embedded mesh vertices. That is, the visual model with detailed geometry does not match the mechanical model on which the modal analysis is performed. The motion of the embedding uses a trilinear interpolation of the mechanical degrees of freedom (DOFs), so we can nevertheless compute the motion of any point on the surface given the mode shapes.

Sound Generation

When rendering the sound with a modal synthesis approach, we do not solve the emission problem, but instead we consider the sound to be simply a sum of damped sinusoids. The activation of this model depends on where the object is hit. If we hit the object at a vibration node of a mode, then that mode will not vibrate, but others will. This is what we refer to as the Sound Map, which could also be called a sound excitation map as it indicates how the different modes are excited when the object is struck at different locations.

The sound resulting from an impact on a specific location on the surface is calculated as a sum of n damped oscillators:

s(t) = 1 n ai sin(wit)e -d i t (1)
where wi, di, and ai are respectively the frequency, the decay rate and the gain of the mode i.

In our method, we synthesize the sounds via a reson filter (see, for example, Van den Doel et al. [START_REF] Kees Van Den | Foley automatic: physically-based sound effects for interactive simulation and animation[END_REF]). This choice is made based on the effectiveness for real-time audio processing. No radiation properties are considered; our study focuses specifically on effective modal synthesis. However, radiation can be computed in a number of ways [START_REF] James | Synthesizing sounds from physically based motion[END_REF][START_REF] James | Precomputed acoustic transfer: output-sensitive, accurate sound generation for geometrically complex vibration sources[END_REF]. As the motions of objects are computed with modal analysis, surfaces can be easily analyzed to determine how the motion will induce acoustic pressure waves in the surrounding medium. Finally, our study does not consider contactposition dependent damping or changes in boundary constraints, as might happen during moments of excitation. Instead we use a uniform damping value for the sounding object.

Position Dependent Sound Rendering

To properly render impact sounds of an object, the method must preserve the sound variety when hitting the surface at different locations. For example, consider the metal bowl, modeled by a triangle mesh with 274 vertices, shown in Figure 1.

Figure 1: A sounding metal bowl: sound synthesis is performed for excitation on specific locations on the surface: points 30, 40 and 52.

Figure 1 specifies where the bowl is hit. We take 3 different locations, i.e., top, side and bottom, on the surface of the object where the object is impacted. The excitation force is modelled as a dirac, such as a regular impact. The material of the bowl is aluminium, with the parameters 69×10 9 for Young Modulus, 0.33 for Poisson coefficient, and 2700 kg/m3 for the volumic mass. The Rayleigh damping parameters for stiffness and mass are set to 3×10 -7 and 10. The use of a constant damping ratio is a simplification that still produces good results.

We compare our approach to modal analysis performed first using tetrahedralization with Tetgen 3 with 822 modes. Our method uses hexahedral finite elements and is applied with a grid of 2×2×2 cells, leading to 81 modes. However, to adapt the stiffness of a cell according to its content, the mesh is refined more precisely than desired for the animation. The information is propagated from fine cells to coarser cells. For this example, the elements of the 2×2×2 coarse grid resolution approximates mechanical properties propagated fom a fine grid of 4×4×4 cells.

The frequency content of the sound resulting from impact at the 3 locations on the surface is shown in Figure 2.

Figure 2: Sound synthesis with a modal approach using classical tetrahedralization with 822 modes (left) and our method with a 2×2×2 hexahedral FEM resolution, leading to 81 modes (right): power spectrum of the sounds emitted when impacting at the 3 different locations shown in Figure 1, (from top to bottom) points 30, 52 and 40.

In Figure 2, each power spectrum is normalized with the maximum amplitude in order to factor out the magnitude of the impact. The eigenvalues that correspond to vibration modes will be nonzero, but for each free body in the system there will be six zero eigenvalues for the body's six rigid-body freedoms. Only the modes whith nonzero eigenvalue are kept. Thus, 816 modes are finally used for sound rendering with the tetrahedralization method and 75 with our hexahedral FEM method.

The movie provided4 compares the sounds synthesized with the tetrahedral FEM and the hexadedral FEM approaches. While Figure 2 highlights the visual differences in the frequency content, we notice in listening to the synthesized sounds that those generated by our method are quite similar to those created with the standard tetrahedralization, even when significantly fewer vibration modes are used (i.e., 75 in contrast to 816).

ROBUSTNESS AND MULTI-SCALE RESULTS

Computing modes for complex geometries can become prohibitively expensive especially when numerous sounding objects have to be processed. As an example, the actual cost of computing the partial eigenvalue decomposition using a tetrahedralization in the case of a bowl with 274 vertices and generating 2426 tetrahedras is 5 minutes with an Intel Core Duo with 2.33 GHz and 2 GB of memory. The number of tetrahedras determine the dimension of the system to solve. To avoid this expense, we provide a method that greatly simplifies the modal parameter extraction even for non-manifold geometries that include both volumetric and surface parts. Our technique consists of using multi-resolution hexahedral embeddings.

Robustness

Most approaches for tetrahedral mesh generation have limitations. In particular, an important requirement imposed by the application of deformable FEM is that tetrahedra must have appropriate shapes, for instance, not too flat or sharp. By far the most popular of the tetrahedral meshing techniques are those utilizing the Delaunay criterion [START_REF] Richard | A condition guaranteeing the existence of higher-dimensional constrained delaunay triangulations[END_REF]. When the Delaunay criterion is not satisfied, modal analysis using standard tetrahedralization is impossible. In comparison with tetrahedralization methods, our technique can handle complex geometries and adequatly performs modal analysis. Figure 3 gives an example of problematic geometry for tetrahedralization because of the presence of very thin parts, specifically the blades that protrude from either side. Figure 3: An example of a complex geometry that can be handled with our method. The thin blade causes problems with traditional tetrahedralization methods.

We suppose the object is made of aluminium (see Section 3.3 for the material parameters). We apply a coarse grid of 7×7×7 cells for modal analysis. The coarse level encloses the mechanical properties of a fine grid of 14×14×14 cells. Figure 5 shows the power spectrum of the sounds resulting from impacts, modelled as a dirac, on 5 different locations. Each power spectrum is normalized with the maximum amplitude of the spectrum in order to factor out the magnitude of the impact.

Figure 5 shows that the Sound Map is preserved; we can observe that the different modes have varying amplitude depending on the location of excitation. It is interesting to examine the quality of the sound rendered when hitting the wings. Because this part is lightweight compared to the rest of the object, the amplitude of higher frequencies is more pronounced than at other locations.

A Multi-scale Approach

To study the influence of the number of hexahedral finite elements on the sound rendering, we model a sounding object with different resolutions of hexahedral finite elements. We have created a squirrel model with 999 vertices which we use as our test sounding object. Its material is pine wood, which has parameters 12×10 9 for Young Modulus, 0.3 for Poisson coefficient, 750 kg/m 3 for The coarse grid resolution for finite elements is set to 2×2×2, 3×3×3, 4×4×4, 8×8×8 and 9×9×9 cells. In this example, the finer grid resolution is one level up to the one of coarse grid, that is, a coarse grid of 2×2×2 cells has a fine level of 4×4×4 cells.

Results show that the frequency content of sounds depend on the location of excitation and on the resolution of the hexahedral finite elements. The higher resolution models have a wider range of frequencies because of the supplementary degrees of freedom. We also observe a frequency shift as the FEM resolution increases. Note that a 2×2×2 grid represents an extremely coarse embedding, and consequently it is not surprising that the fraquency content is different at higher resolution. Nevertheless, there are still some strong similarities at the dominant frequencies. Above all, an important feature is the convergence in frequency content as the FEM increases. According to Figure 6, a grid of 4×4×4 cells may be sufficient to properly render the sound quality of the object.

Limitations

The Sound Map is influenced by the resolution of the hexahedral finite elements. This is related to the way stiffnesses and masses of different elements are altered based on their contents. As a consequence, a 2×2×2 hexahedral FEM resolution would show much less expressive variation than higher FEM resolution. This is shown in the movie 5 . One approach to improving this would be to use better approximations of the mass and stiffness of coarse elements [START_REF] Nesme | Preserving topology and elasticity for embedded deformable models[END_REF].

Nevertheless, based on the quality of the resulting sounds, and given that increased resolution for the finite elements implies higher memory and computational requirements for modal data, finite elements resolution can be adapted to the number of sounding objects in the virtual scene.

DISCUSSION

Table 1 gives the computation time and the memory usage of the modal data when computing the modal analysis with different FEM resolution on the squirrel model. In this example, the finer grid resolutin is one level up to the one of coarse grid, that is, a coarse grid of 4×4×4 cells has a fine level of 8×8×8 cells. putation times of our unoptimized initial implementation on a 2.33 GHz Intel Core Duo. Despite the fact that audio is considered a very important aspect in virtual environments, it is still considered to be of lower importance than graphics. We believe that physically modeled audio brings a significant added value in terms of realism and the increased sense of immersion.

The use of physics engines is becoming much more widespread for animated interactive virtual environments; the interface between these engines and audio has often been one of the obstacles for the adoption of physically based sound synthesis in simulations. This is often due to the lack of appropriate design choices in the two interfaces that prevent them from working together effectively.

Our method is built on a physically based animation engine, Sofa Framework. As a consequence, problems of coherence between physics simulation and audio are avoided by using exactly the same model for simulation and sound modeling. 5 Additional material: http://www-sop.inria.fr/reves/Cecile.Picard

CONCLUSION

We propose a new approach to modal analysis using automatic voxelization of a surface model and automatic tuning of the finite elements parameters, based on the distribution of material in each cell. Our goal is to perform sound rendering in the context of an animated real-time virtual environment, which has specific requirements, such as real-time processing and efficient memory usage.

We have shown that in simple cases our method gives similar results as traditional modal analysis with tetrahedralization for simple cases. For more complex cases, our approach provides plausible results. In particular, sound variety along the object surface, the Sound Map, is well preserved.

Our technique can handle complex non-manifold geometries that include both volumetric and surface parts, which can not be handled by previous techniques. We are thus able to compute the audio response of numerous and diverse sounding objects, such as those used games, training simulations, and other interactive virtual environment.

Our solution allows a multi-scale solution because the number of hexahedral finite elements only loosely depends on the geometry of the sounding object.

Finally, since our method is built on a physics animation engine, the Sofa Framework, problems of coherence between simulation and audio can be easily addressed, which is of great interest in the context of interactive environment.

In addition, due to the fast computation time, we are hopeful that real-time modal analysis will soon be possible on the fly, with sound results that are approximate but still realistic for virtual environments.

Figure 4 :

 4 Figure 4: Test impacts for sound generation are simulated on 5 different locations on the surface of the complex geometry: points 100, 135, 146, 200 and 400.

Figure 5 :

 5 Figure 5: The power spectrum of the sounds resulting from impacts at the 5 different locations shown in Figure 4: (from top to bottom) points 100, 135, 146, 200 and 400. Note that the audible response is different based on where the object is hit.

Table 1 :

 1 These are com-Computation time and memory usage for different grid resolutions.

	Coarse Grid Resolution Computation Time Memory Usage
	(cells)	(seconds)	(MB)
	7×7×7	115.11	11.309
	6×6×6	39.14	5.698
	5×5×5	12.98	2.663
	4×4×4	5.35	1.042

Simulation Open Framework Architecture; http://www.sofa-framework.org/ SOFA scenes are modeled using a data structure similar to hierarchical scene graphs which allows the physical objects to be split easily into collections of independent components, each describing one feature of the model. Moreover, simulation algorithms are also modeled as components in the scene graph, providing us with the same flexibility for algorithms as for models.Elastic deformations are used to generate the audio signal. Before performing modal decomposition, we must first select a deformable modeling method that can be used to generate the stiffness and the mass matrices of the mechanical system. A variety of methods could be used, including particle systems[START_REF] Raghuvanshi | Interactive sound synthesis for large scale environments[END_REF] or finite differences methods. The tetrahedral finite element method has also been used[START_REF] James | Synthesizing sounds from rigid-body simulations[END_REF]. However, tetrahedral meshes are computationally expensive for complex geometries, and can be difficult to tune. As an example, in the tetrahedral mesh

generator Tetgen 2 , the mesh element quality criterion is based on the minimum radius-edge ratio, which limits the ratio between the radius of the circumsphere of the tetrahedron and the shortest edge length.Our method is inspired from work by Nesme et al.[START_REF] Nesme | Animating shapes at arbitrary resolution with non-uniform stiffness[END_REF]. It uses hexahedral finite element for computing the mass and stiffness matrices of the mechanical system. The technique can be summarized as follows. An automatic high-resolution voxelization of the geometric object is first built. The voxelization initially concerns the surface of the geometric model, while the interior is automatically filled when the geometry represents a solid object. The voxels are then recursively merged up to an arbitrary coarser mechanical resolution. The merged voxels are used as hexahedral finite elements embedding the detailed geometrical shape. At each level, the mass and stiffness of a merged voxel are deduced from its eight children, using a weigthed average that takes into account the distribution of material. With this method, we can handle objects with geometries that simultaneously include volumetric and surface parts; thin or flat features will occupy voxels and will thus result in the creation of mechanical elements that approximate their shape (see Section 4.1).We extend the method for microscopic deformations that allows sound rendering. Thus, in order to compute the modal parameters, we compute the assembled mass and the assembled stiffness matrices for the object by summing the contribution of each cell. Then, we solve the decoupled system to extract the modal parameters as explained in Section 2. Our preprocessing step that performs modal analysis can be summarized as follows.

http://tetgen.berlios.de/

Tetrahedral Mesh Generator: http://tetgen.berlios.de/

Additional material: http://www-sop.inria.fr/reves/Cecile.Picard

Video games studio; http://www.eden-games.com/

Proc. of the 12 th Int. Conference on Digital Audio Effects (DAFx-09), Como, Italy, September[START_REF] James | Synthesizing sounds from rigid-body simulations[END_REF][START_REF] Adrien | The missing link: modal synthesis[END_REF][START_REF] Perry | Real Sound Synthesis for Interactive Applications[END_REF][START_REF] Iovino | Recent work around modalys and modal synthesis[END_REF] 2009

ACKNOWLEDGMENTS

This work was partly funded by Eden Games 6 , an ATARI Game Studio in Lyon, France. We would like to thank Nicolas Tsingos for his input on an early draft.

Pt 29

Pt 287 Pt 986