Juliette Chabassier

Marc Duruflé

P Joly

TIME DOMAIN SIMULATION OF A PIANO. PART 2 : NUMERICAL ASPECTS

Keywords: 1991 Mathematics Subject Classification. 00A71, 00A65, 65P05, 65N25, 35Q72, 35L05 piano model, energy preserving schemes, numerical methods

HAL is a

Introduction

This article is the second of a series of two papers devoted to the numerical simulation of the piano. In the first part [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF], we have constructed our piano model as a coupled system of PDEs, based on a complete description of the physical mechanisms involved in the generation of a piano sound from the hammer blow to the sound radiation in the air, through the string vibrations, and their transmission to a radiating soundboard. This second part is devoted to the construction and analysis of a numerical method for the approximation of this system, and to the presentation of numerical results obtained with the corresponding computation code. Following the pioneering work by [START_REF] Giordano | Physical modeling of the piano[END_REF], this work is a step forward in the modeling of the piano, regarding modeling issues, as it has been shown in [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF], as well as numerical methods, as it will be demonstrated in this paper. Authors have explored many possibilities offered by numerical methods in order to better understand the physics of the different parts of the piano. For instance, numerical methods for the soundboard can be found in [KW87, Gio97, BIJ03, OBCQ08, EBR13] and the down bearing is simulated in [START_REF] Mamou-Mani | Numerical simulation of a piano soundboard under downbearing[END_REF], numerical methods for the string and eventually its coupling with the hammer head are explored in [Bou88, CA94a, CA94b, BBKMI03, BS04, Bil05], numerical methods for the flexible hammer shank are available in [START_REF] Izadbakhsh | Dynamic modeling and experimental testing of a piano action mechanism with a flexible hammer shank[END_REF]. In this context, our approach in this second paper aims at designing a comprehensive and innovative numerical tool which considers the piano as a whole physical system.

Although the reading of article [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF] would be very helpful for the reader's understanding, we have chosen to review briefly in section 1 the equations of the piano model. In section 2 we present the variational formulation of the problem, which is a basis for its space discretisation presented in section 3. The main section of this work is section 4, devoted to time discretisation, which contains most original and tricky aspects of our method.

Due to the complexity of the global model, we have chosen to present separately the schemes used for each physical subsystems. The most delicate step concerns the discretisation of the nonlinear strings equations, and is presented in section 4.2. The hammer strings coupling is treated in section 4.3. The strings soundboard coupling is dealt with in section 4.4 and the approximation of the soundboard acoustics interaction is described in section 4.5. At last, in section 4.6, we present our numerical scheme for the global piano model. Section 5 presents the practical and algorithmic aspects of the work that are used for the implementation of our computational code. Finally, we present various numerical results that illustrate the ability of our computational code to describe the full behavior of the piano (section 6.1), show that our piano model does fulfill the objectives previously assigned in [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF] (section 6.2) and demonstrate that we can reproduce with very good qualitative property the experimental property (section 6.3). This is to our knowledge the first comprehensive simulations of the piano that achieve this degree of realism.

In our presentation, we shall insist on the stability issues related to the method, which constitutes one of the main difficulties of the numerical analysis. This will be systematically treated via energy methods. Even more, we shall show how the desire to ensure good energy properties for our numerical method has guided the construction of the schemes.

Review of the piano model

In this first section, we will recall the model that we have retained for the piano. This summarizes the content of the article [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF] to which we refer the reader for additional details and information. However, we try to give the necessary amount of information to keep this article self-contained. Also, for facilitating the presentation of our numerical method, we found useful to introduce some abstract notation that helps for a more compact writing of the equations.

In what follows, we shall use the 3D space coordinates (x, z), x = (x 1 , x 2 ) and denote (e 1 , e 2 , e z ) the corresponding orthonormal basis. All our notations are illustrated in figure 1. The piano, strings apart, is represented as the 3D object Ω f ∪ ω where

• ω is the 2D domain in the plane z = 0 that represents the soundboard, seen as a thin plate,

• Ω f (where Ω f is a bounded open set of R 3 ) represents the rim of the piano, which will be considered as an obstacle to the sound propagation.

The 3D domain occupied by the air is the exterior domain Ω (B is a box in R 3 containing Ω f ∩ ω with outward normal n B )

Ω = B \ Ω f ∩ ω . (1)
We shall distinguish, for the same geometric object ω, the two sides of the soundboard by introducing

ω + ≡ ∂Ω + ∩ ω, ω -≡ ∂Ω -∩ ω, ,
where

Ω ± = Ω ∩ {± z > 0} (2)
The boundary of Ω is therefore ∂Ω = ∂Ω f ∪ ω + ∪ ω -and a function q defined in Ω (in H 1 (Ω)) will have two distinct traces on ω, namely q| ω + and q| ω -.

Let us simply mention that we are interested in the coupling of a chord of strings : several parallel strings are hit by the same hammer, N s is the number of strings (1,2 or 3 in practice) i is an index for the strings. We denote x the abscissa along the strings, with length L. We assume that the motion of the strings is confined to a plane orthogonal to the soundboard.

The equations of the piano model are given by the following system.

                                                                                                                                           Find ξ : R + -→ R, q i : [0, L] × R + -→ R N , 1 ≤ i ≤ N s , U p : ω × R + -→ R P , N = 3, P = 3, U a = (p, V ) : Ω × R + -→ R × R 3 such that Hammer's equations          m H d 2 ξ dt 2 (t) = - i F i (t), F i (t) = k H Φ H e i (t) + r H d dt Φ H e i (t) , e i (t) = ξ -ξ i (t) + ξ(t) + , ξ i (t) = L 0 q i • ν s (x, t) δ H (x -x i ) dx, (3a)
Strings equations

M ∂ 2 t q i + ∂ t R q i -∂ x (Γ ∂ x q i ) -∂ x (∇ p H i (∂ x q i , q i )) + ∇ q H i (∂ x q i , q i ) = F i (t) δ H (x -x i ) ν s (3b)
Strings' boundary conditions and strings / soundboard coupling equations

       q i (L, t) • τ s α = 0, q i (L, t) • ν s α = ω χ ω (x -x a ) U p (x, t) • ν p Jq i (0, t) = 0, Id -J Γ ∂ 2 xt q i + ∇ p H i (∂ x q i , q i ) (x e , t) = 0, x e = 0 or L (3c)
Soundboard's equations

M p ∂ 2 t U p + f d (A p ) ∂ t U p + A p (x, ∇ x )U p = χ ω (x -x a ) i ∇ p H i (∂ x q i , q i ) + Γ ∂ 2 xt q i • τ s α + p ω ν p , (3d)
Soundboard's boundary conditions B p (n, ∇ x )U p = 0, on ∂ω, t > 0.

Sound propagation

ρ a ∂ t V + ∇p = 0, µ a ∂ t p + div V = 0, (3f)
Acoustic boundary conditions -soundboard / air coupling conditions

V • e z | ω -= V • e z | ω + = ∂ t U p • ν p , on ω, V • n f = 0 on ∂Ω f and V • n B = 0 on ∂B. (3g)
These equations are completed by initial conditions which simply state that the piano (and the surrounding air) is at rest at t = 0, the instant where the hammer hits the strings with initial velocity v H :

                     ξ(0) = -ξ, dξ dt (0) = v H , q i (x, 0) = ∂ t q i (x, 0) = 0, x ∈ [0, L],
U p (x, 0) = ∂ t U p (x, 0) = 0, x ∈ ω, p(x, z, 0) = 0, V (x, z, 0) = 0, (x, z) ∈ Ω.

(4)

ω - ω + Ω x 1 x 2 z Ω f
Φ H (d) = d ν , for some ν > 1 (5)
• q i (x, t) is the vector of unknowns associated with the string number i and R N is the space for string's unknowns. More precisely, q i (x, t) = (u i , v i , ϕ i )(x, t), where u i is the transverse displacement, v i is the longitudinal displacement and ϕ i is the angle of rotation of the cross sections. The vectors τ s α , ν s α and ν s are the vectors of R N defined by

τ s α = (-sin α, cos α, 0) t , ν s α = (cos α, sin α, 0) t , ν s = (1, 0, 0) t (6)
where α is the angle between the strings and the plane {z = 0}. The function δ H is a regularized delta function and x i is the abscissa of the contact point between the hammer and the i th string. The matrix J is given by diag(1, 1, 0). The N ×N matrices R and Γ are diagonal and positive, M = diag(ρ A, ρ A, ρ I) and the function H i (p, q) : R N × R N → R + is given by:

H i (p, q) = H i,2 (p, q) + U i (p), where (7)                          H i,2 (p, q) = 1 2 A i p • p + C q • q + 2 B q • p (8a)
A i =   T i 0 + AGκ 0 0 0 EA 0 0 0 EI   , B =   0 0 -AGκ 0 0 0 0 0 0   , C =   0 0 0 0 0 0 0 0 AGκ   , (8b)
U i (p) = (EA -T i 0 ) p 2 1 2 + (1 + p 2 ) -p 2 1 + (1 + p 2 ) 2 , p = (p 1 , p 2 , p 3 ), (8c)
where T i 0 is the tension of the string at rest, A and I are respectively the area and inertia momentum of the string's cross section, E and G are the Young and shear moduli of the material, ρ is the density, and κ is the Timoshenko's parameter (see [START_REF] Hutchinson | Shear coefficients for timoshenko beam theory[END_REF] for a discussion on its physical meaning and one of its possible mathematical derivations). Note that we have assumed (this is not essential anyway) that, in a given choir of strings, only the tension depends on the string (which corresponds physically to a slight detuning between them).

• U p (x, t), x ∈ ω is the vector of plate's unknowns associated with the soundboard and R P is the space for plate's unknowns. More precisely, U p = (u p , θ p ) , where u p is the transverse displacement of the plate, and where θ p = (θ 1,p , θ 2,p ) are the two deflection angles representing the rotations of the normal fibers of the plate. Moreover the diagonal matrix M p is given by

M p (x) = diag ρ p (x) δ(x), ρ p (x) δ(x) 3 /12, ρ p (x) δ(x) 3 /12 (9)
where ρ p and δ are the variable density and thickness of the plate. The differential operator A p (x, ∇ x ) is associated with the Reissner-Mindlin system:

A p (x, ∇ x )   u p θ p   =    -div δ(x) κ(x) 2 G(x) ∇ x u p + θ p ) -Div δ 3 (x) 12 C(x) ε(θ p ) + δ(x) κ 2 (x) G(x) ∇ x u p + θ p    (10)
where κ(x) is the shear correction factor of the Reissner-Mindlin model, div (respectively Div) is the 2D divergence of a vector field (respectively a tensor field), ε is the usual 2D deformation tensor, G(x) and C(x) are the usual stiffness tensors of the Reissner-Mindlin model : G(x) is a positive symmetric semi definite operator in R 2 and C(x) is a positive definite operator in S 2 , the set of real 2×2 symmetric matrices. Moreover

B p (n, ∇ x )U p = B p (n, ∇ x )   u p θ p   =   u p C(x) ε(θ p )n   ( 11
)
where n is the normal to ∂ω in the plane z = 0.

Given a damping function f d : R + → R + of the form f d (λ) = αλ + β √ λ + γ, the operator f (A p ) refers to the selfadjoint operator A p associated with the differential operator A p (x, ∇ x ) with boundary conditions B p (n, ∇ x )U p = 0 in the Hilbert space H p = L 2 (ω) P , equipped with the scalar product

U p , U p Hp = M p (x) U p (x) • U p (x) dx. ( 12
)
Introducing a basis of eigenfunctions and associated positive eigenvalues (λ 1 ≤ λ 2 ≤ ...)

   A p (x, ∇ x )W n = λ n M p (x) W n , x ∈ ω, B p (n, ∇ x )W n = 0, x ∈ ∂ω, (13)
chosen to be orthonormal in H p , the operator f d (A p ) is the self adjoint operator diagonalized in this basis so that

∀ n ∈ N * , f d A p W n = f d (λ n )W n . (14)
The function χ ω (x) is again a 2D regularized delta function, and x a is the position of the attach point of the chord of strings on the soundboard. The vector ν p ∈ R P is equal to (1, 0, 0) t . The jump of pressure [p] ω across ω is defined by p| ω --p| ω + .

• p(x, z, t), (x, z) ∈ Ω is the acoustic pressure and V (x, z, t), (x, z) ∈ Ω is the acoustic velocity field. The coefficients ρ a and µ a represent the density and Lamé's coefficient of the air.

Variational formulation

In this section, we present the weak formulation, or variational formulation in space, that will be the basis of the space discretisation of the problem. The first two boundary conditions in (3c) are essential conditions that could be included in the variational spaces for the string's and plate's unknowns. However, especially for the second equation, this introduces a strong coupling between the two sets of unknowns, which is a priori incompatible with the choice of discretizing plate and string unknowns in a very different way, in both space and time. That is why, in order to ensure a better decoupling, we have chosen to treat these conditions in a weak form by introducing two scalar unknowns per string that can be interpreted as Lagrange multipliers associated with these constraints (physically, F p i is the vertical component of the force exerted by the ith string on the plate, while F r i is the horizontal one):

F p i = ∇ p H i (∂ x q i , q i ) + Γ ∂ 2 xt q i • ν s α (15a) F r i = ∇ p H i (∂ x q i , q i ) + Γ ∂ 2 xt q i • τ s α (15b)
In the formulation below, most unknowns are seen as functions of time with values in appropriate Hilbert spaces. That is why we introduce the following functional spaces, respectively for the string, plate and acoustic unknowns :

V s = {q = (u, v, ϕ) ∈ H 1 (0, L) 3 / u(0) = v(0) = 0 ⇐⇒ Jq(0) = 0} ⊂ H s = L 2 (0, L) 3 V p = {U p = (u p , θ) ∈ H 1 (ω) 3 / u p | ∂ω = 0 ⇐⇒ U p • ν p | ∂ω = 0} ⊂ H p = L 2 (ω) 3 , V a = H 1 (Ω) × L 2 (Ω) 3 ⊂ H a = L 2 (Ω) × L 2 (Ω) 3 . ( 16
)
The space V is defined as the following product

V = (V s ) Ns × V p × V a (17)
In the space V s , we define the symmetric positive bilinear forms

m s (q i , q i ) = L 0 M q i • q i dx, a s d (q i , q i ) = L 0 R q i • q i dx + L 0 Γ ∂ x q i • ∂ x q i dx (18)
as well as the semi-linear (linear in the second argument, nonlinear in the first one) form

a s i (q i , q i ) = L 0 ∇ p H i (∂ x q i , q i ) • ∂ x q i dx + L 0 ∇ q H i (∂ x q i , q i ) • q i dx (19)
According to the decomposition (7), we note that a s i (q i , q i ) can be split into two contributions, which permits to identify the linear part and nonlinear part of the equation :

a s i (q i , q i ) = a s,2 i (q i , q i ) + a s,nl i (q i , q i ) (20)
where

         a s,2 i (q i , q i ) = L 0 ∇ p H i,2 (∂ x q i , q i ) • ∂ x q i dx + L 0 ∇ q H i,2 (∂ x q i , q i ) • q i dx (21a) a s,nl i (q i , q i ) = L 0 ∇U i (∂ x q i ) • ∂ x q i dx. ( 21b
)
Finally we introduce the linear forms

L H ( q i ) = L 0 δ H (x -x i ) q i • ν s dx, L p ( q i ) = q i (L) • ν s α , L r ( q i ) = q i (L) • τ s α . (22)
In the space V p , we define the symmetric positive bilinear forms

         m p (U p , U p ) = L 0 M p U p • U p dx, a p (U p , U p ) = ω δ 3 (x) 12 C(x) ε(θ p ) : ε( θ p ) + ω δ(x) κ 2 (x) G(x) ∇u p + θ p • ∇ u p + θ p dx (23)
where ε : ε denotes the usual scalar product in S 2 , the space of 2 × 2 symmetric matrices.

a p d (U p , U p ) = m f (λ m ) m p (U p , W m ) m p ( U p , W m ) (24)
Finally, we define

L s ( U p ) = ω χ ω (x -x a ) U p • ν p dx. (25)
In the space V a , we define the bilinear forms (respectively symmetric positive definite and skew symmetric)

m a (U a , U a ) = Ω µ a p p + ρ a V • V dx, b a (U a , U a ) = Ω ∇p • V -∇ p • V dx, (26)
Finally, we define the bilinear form on V p × V a (we recall that U a = ( p, V ))

c ap ( U p , U a ) = ω [ p ] ω U p • ν p dx (27)
The variational formulation of system (3) is

                                                               Find ξ, (F p i ) 1≤i≤Ns , (F r i ) 1≤i≤Ns , (q i ) 1≤i≤Ns , U p , U a ) : R + -→ R × R Ns × R Ns × V such that for any ( q i ) 1≤i≤Ns , U p , U a ) ∈ V      m H d 2 ξ dt 2 (t) = - i F i (t), F i (t) = k H Φ H e i (t) + r H d dt Φ H e i (t) , e i (t) = ξ + ξ i (t) -ξ(t) + , ξ i (t) = L H (q i ) (28a) d 2 dt 2 m s (q i , q i ) + d dt a s d (q i , q i ) + a s i (q i , q i ) = F i (t) L H ( q i ) + F p i (t) L p ( q i ) + F r i (t) L r ( q i ) (28b) L p (q i ) = L s (U p ), L r (q i ) = 0 (28c) d 2 dt 2 m p (U p , U p ) + d dt a p d (U p , U p ) + a p (U p , U p ) = - i F p i (t) L s ( U p ) + c ap ( U p , U a ) (28d)
d dt m a (U a , U a ) + b a (U a , U a ) = - d dt c ap (U p , U a ) (28e)
The derivation of (28) from (3) is pretty standard. Let us simply make somme comments.

• Equation (28b) is obtained from (3b) after multiplication (inner product in R N ) by a test field q i ∈ V s and integration between 0 and L. After integration by parts, most of the boundary terms vanish because J q i (0) = 0 ( q i ∈ V s ) and thanks to (3c) (second line). The remaining term is obtained from the definition of F p i and F r i (see (15)) and the definition of L p and L r (see ( 22)). • Equation (28d) is obtained from (3d), after multiplication (inner product in R M ) by a test field U p ∈ V p and integration (by parts) over ω.

• To obtain (28e) we first multiply (inner product in R 3 ) the first equation by a test field V and integrate over Ω. Next we multiply the second equation of (3f) by a test function p and integrate over Ω. Adding the two equalities and using Green's formula we obtain (with obvious notation, n f denotes the unit vector normal to ∂Ω, outward with respect to Ω)

d dt m a (U a , U a ) + b a (U a , U a ) = ∂Ω f p V • n f + ω + p V • e z - ω - p V • e z + ∂Ω B p V • n B
One concludes using (3g).

It is easy to recover formally the energy identity satisfied by any smooth enough solution of (28) by choosing

q i = dq i dt , U p = dU p dt , U a = U a .
This leads to

d dt E tot (ξ, q i , U p , U a ) = - i L 0 R ∂ t q i • ∂ t q i - i L 0 Γ ∂ 2 x,t q i • ∂ 2 x,t q i -a p d (∂ t U p , ∂ t U p ) - i r H Φ ′ H e i ∂ t e i 2 ≤ 0 (29)
In (29), with Ψ H (d) = d 0 Φ H (s) ds ≥ 0, the total energy of the system is given by

E tot (ξ, q i , U p , U a ) = m H 2 ∂ t ξ 2 + i k H Ψ H e i + i E s (q i ) + E p (U p ) + E a (U a ). ( 30
)
where E s (q i ) is the energy of the i th string, E p (U p ) is the soundboard energy and E a (U a ) the acoustic energy:

               E s (q i ) = 1 2 m s (∂ t q i , ∂ t q i ) + L 0 H i (∂ x q i , q i ) dx (31a) E p (U p ) = 1 2 m p (∂ t U p , ∂ t U p ) + 1 2 a p (U p , U p ) (31b) E a (U a ) = 1 2 m a (U a , U a ) (31c)
The only non trivial point in the derivation of (29) is

a s i (q i , ∂ t q i ) = L 0 ∇ p H i (∂ x q i , q i ) • ∂ 2 x,t q i + ∇ q H i (∂ x q i , q i ) • ∂ t q i dx = d dt L 0 H i (∂ x q i , q i ) dx
as a consequence of the chain rule.

Space discretisation and algebraic formulation of the semi-discrete problem

Construction of the semi-discrete problem by a Galerkin method

The space discretisation of the variational problem (28) will be based on

• one dimensional higher order Lagrange finite elements for the approximation of the strings' unknowns q i ,

• hexahedral continuous Lagrange higher order finite elements for the pressure field p and discontinuous Lagrange higher order finite elements for the acoustic velocity V ,

• a truncated modal (or spectral) approximation for the plate unknowns.

We introduce as an abstract approximation parameter devoted to tend to zero:

h = ∆x, M -1 , h (32)
where

• ∆x is the step size of a uniform grid of [0, L],

• M is an integer that represents a spectral truncation order,

• h represents the mesh size of a three dimensional mesh of Ω.

More precisely, for the strings, given a polynomial degree r s , we consider the finite dimensional subspace of V s :

V s,∆x = q i,h ∈ V s such that ∀ ℓ < L/∆x, q i,h [ℓ∆x,(ℓ+1)∆x] ∈ P rs ⊂ V s . ( 33
)
where

P rs = Span{x i , 1 ≤ i ≤ r s }.
For the plate unknowns we define

V p,M = Span{W m , 1 ≤ m ≤ M } ⊂ V p ( 34
)
where W m are the plate eigenmodes defined by (13). For the acoustic unknowns, one uses a quasi uniform hexahedral mesh of Ω with step size h > 0 (of course, this implies in practice a geometrical approximation of the domain Ω but, for simplicity, we shall omit these details here) :

T h (Ω) = {K i } 1≤i≤Na , Ω = Na i=1 K i and K i ∩ K j = ∅, if i = j, h = max i diam(K i ) (35)
where K i is a hexahedron defined as the image of the unit cube by a trilinear map Φ i . Given a polynomial r a > 0, we consider the discrete space

V a,h = U a,h = (p h , V h ) ∈ V a , such that p h Ki • Φ i ∈ Q ra and V h Ki • Φ i ∈ Q 3 ra ⊂ V a (36)
where

Q ra = Span{x i 1 x j 2 z k , 1 ≤ i, j, k ≤ r a }.
This suggests to define the discrete approximation space for V (see (17)):

V h = (V s,∆x ) Ns × V p,M × V a,h (37)
and the natural semi-discrete problem consists in rewriting the variational problem (28), simply by replacing V by V h . However, in practice, we are led to an additional approximation step due to the following reasons (i) We use systematically quadrature formula (element per element) for evaluating all integrals, over [0, L] or Ω, appearing in the formulation. More precisely, for the 1D integrals along each segment of the string, we use the unique Gauss-Legendre formula with r s + 1 points that integrates exactly polynomials of degree 2 r s + 1. This presents the interest of preserving the accuracy provided by the interpolation properties of V s by V s,∆x . For the volume integrals in Ω, all integrals on hexahedra are evaluated by using a quadrature formula on the unit cube obtained by a tensor product of 1D Gauss-Lobatto formula of order r a in each variable. By using Gauss-Lobatto points both for integration and interpolation, mass lumping is achieved, i.e. the mass matrix associated with the bilinear form m a is diagonal. We refer the reader to [GC00] for more details .

(ii) The eigenmodes W m are not known analytically and must be computed approximately via numerical methods. We have chosen to use higher order finite elements methods based on a quadrilateral mesh T h (ω) = {K ω,ℓ , 1 ≤ ℓ ≤ N p } of the domain ω which coincides for simplicity with the intersection between the 3D mesh T h (Ω) and ω. Each quadrilateral K ω,ℓ is the image of the unit square by a bilinear transformation Φ ℓ . Given a polynomial degree r p , we define

V p,h = U p,h ∈ V p , such that U p,h K ℓ,ω • Φ ℓ ∈ Q P rp ( 38
)
where

Q r = Span{x i 1 x j 2 1 ≤ i, j ≤ r p }.
We then solve the eigenvalue problem

Find W h ∈ V p,h \ {0} and λ h ∈ R + , such that for all U p,h ∈ V p,h , a p h (W h , U p,h ) = λ h m p h (W h , U p,h ) (39)
The notations a p h and m p h refer to the use of Gauss-Lobatto quadrature formulas which are the 2D equivalent of the 3D formula described in point (i). We order increasingly the positive eigenvalues, and, assuming that dim

(V p,h ) ≥ M , we select the M first eigenvalues {λ h,m , 1 ≤ m ≤ M }. Denoting {W h,m , 1 ≤ m ≤ M } the corresponding orthonormal family of eigenvectors, which are approximations of {W m , 1 ≤ m ≤ M },
we replace the space V p,M with the new space

V p,M,h = Span{W h,m , 1 ≤ m ≤ M } which means that V h = (V s,∆x ) Ns × V p,M,h × V a,h of dimension N s h × M × N a h . (40)
Finally, we approximate the bilinear form

a p d (U p,h , U p,M,h ), for any (U p,h , U p,h ) ∈ V 2 p,M,h , with a p d,h (U p,h , U p,h ) = M m=1 f (λ h,m ) m p h (U p,h , W h,m ) m p h ( U p,h W h,m ). (41)
In the following, as in (39), we shall use systematically the subscript h to indicate the use of quadrature formulas for numerical integration. For instance m s h (p, q) replaces m s (p, q). The semi-discrete problem reads

                                                                                         Find ξ h , (F p i,h ) 1≤i≤Ns , (F r i,h ) 1≤i≤Ns , (q i,h ) 1≤i≤Ns , U p,h , U a,h ) : R + -→ R × R Ns × R Ns × V h such that for any ( q i,h ) 1≤i≤Ns , U p,h , U a,h ) ∈ V h      m H d 2 ξ h dt 2 (t) = - i F i,h (t), F i,h (t) = k H Φ H e i,h (t) + r H d dt Φ H e i,h (t) , e i,h (t) = ξ + ξ i,h (t) -ξ h (t) + , ξ i,h (t) = L H,h (q i,h ) (42a)
d 2 dt 2 m s h (q i,h , q i,h ) + d dt a s d,h (q i,h , q i,h ) + a s i,h (q i,h , q i,h ) = F i,h (t) L H,h ( q i,h ) + F p i,h (t) L p,h ( q i,h ) + F r i,h (t) L r,h ( q i,h ) (42b) L p,h (q i,h ) = L s,h (U p,h ), L r,h (q i,h ) = 0 (42c) d 2 dt 2 m p h (U p,h , U p,h ) + d dt a p d,h (U p,h , U p,h ) + a p h (U p,h , U p,h ) = - i F p i,h (t) L s,h ( U p,h ) + c ap,h ( U p,h , U a,h ) (42d)
d dt m a h (U a,h , U a,h ) + b a h (U a,h , U a,h ) = - d dt c ap,h (U p,h , U a,h ) (42e)

Algebraic form of the semi-discrete problem

Concerning the string's unknowns, we denote Q h and Q h in R N s h be the vector of coordinates of q h and q h ∈ V s,∆x in the Gauss-Lobatto Lagrange finite element basis. For the linear part of the string' s equation, using Riesz's representation theorem, we define the following N s h ×N s h symmetric and positive matrices

M s h , R s h , A s i,h
such that, for any

(q h , q h ) ∈ V 2 s,∆x M s h Q h • Q h = m s h (q h , q h ), A s i,h Q h • Q h = a s i,h (q h , q h ), R s h Q h • Q h = a s d,h (q h , q h ). (43)
In the same way, for the nonlinear part of the equation, we define the application ∇U i : R N s h → R N s h by:

∇U i (Q h ) • Q h = a s,nl i,h (q h , q h ), (44)
Finally, we define the vectors L H,h , L p,h and L r,h in R N s h such that :

L H,h • Q h = L H,h ( q h ) , L p,h • Q h = L p,h ( q h ) , L r,h • Q h = L r,h ( q h ) (45)
For the plate unknowns, if U p,h and U p,h are the vectors of coordinates in the modal basis

{W h,m , 1 ≤ m ≤ M } of U p,h and U p,h ∈ V 2 p,M,h , we observe that m p h (U p,h , U p,h ) = U p,h • U p,h , a p h (U p,h , U p,h ) = Λ p h U p,h • U p,h , a p d,h (U p,h , U p,h ) = Λ p d,h U p,h • U p,h (46)
where the M × M diagonal matrices Λ p h and Λ p d,h are defined by

Λ p h = diag {λ h,m } 1≤m≤M , Λ p d,h = diag {f (λ h,m )} 1≤m≤M , (47)
Moreover, the vector L s,h is defined by

L s,h • U p,h = L s,h ( U p,h ), ∀ U p,h ∈ V p,M,h (48)
For the acoustic unknowns, we denote U a,h and U a,h the vectors of coordinates of U a,h and U a,h ∈ V a,h in the Gauss-Lobatto Lagrange finite element basis. We define the symmetric and positive

N a h × N a h matrices M a h , B a
h and C ap,h such that for any U a,h ∈ V a,h and any

U p,h ∈ V p,M,h , M a h U a,h • U a,h = m a h (U a,h , U a,h ), B a h U a,h • U a,h = b a h (U a,h , U a,h ), (49)
C ap,h U p,h • U a,h = c ap,h U p,h , U a,h , (50)
Remark 1. According to the expression of the bilinear forms (26), note that the matrices M a h and B a h have a particular block decomposition, that will be exploited for time discretization (see section 4.5.1) :

M a h = M p h 0 0 M V h , B a h = 0 C h -C t h 0 ( 51
)
where the blocks correspond to the decomposition of U a,h as U a,h = (p h , V h ). (cf (36)) and correspondingly U a,h = (P h , V h ) t . In the same way, we note that, as the bilinear form c ap (•, •), and thus c ap,h (•, •), only involves the acoustic pressure and not he velocity field (see 27), the coupling matrix C ap,h has the following block structure:

C ap,h = C pp,h 0 (52)
We shall use this decompositions for the time discretization (see section 4.5).

Then, it is straightforward to check that (42) is equivalent to the following nonlinear ordinary differential system

                                                         Find ξ h , (F p i,h ) 1≤i≤Ns , (F r i,h ) 1≤i≤Ns : R + → R × R Ns × R Ns and (Q i,h ) 1≤i≤Ns , U p,h , U a,h ) : R + → R N s h × R M × R N a h such that (53a)      m H d 2 ξ h dt 2 (t) = - i F i,h (t), F i,h (t) = k H Φ H e i,h (t) + r H d dt Φ H e i,h (t) , e i,h (t) = ξ + ξ i,h (t) -ξ h (t) + , ξ i,h (t) = L H,h • Q i,h (53b)
M s h d 2 dt 2 Q i,h + R s h d dt Q i,h + A s i,h Q i,h + ∇U i (Q i,h ) = F i,h (t) L H,h + F p i,h (t) L p,h + F r i,h (t) L r,h (53c)
L p,h • Q i,h = L s,h • U p,h , L r,h • Q i,h = 0 (53d) d 2 dt 2 U p,h + Λ p d,h d dt U p,h + Λ p h U p,h = - i F p i,h (t) L s,h + t C ap,h U a,h (53e)
M a h d dt U a,h + B a h U a,h = - d dt C ap,h U p,h (53f)
completed by initial conditions which are the same as (4) but written for the semi discrete unknowns.

Well-posedness and stability of the semi-discrete problem.

Following the same approach as in the end of section 2, it is immediate to recover formally the energy identity satisfied by any smooth enough solution of (42) by choosing

q i,δx = dq i,h dt , U p,h = dU p,h dt , U a,h = U a,h .
This leads to

d dt E tot (ξ h , q i,h , U p,h , U a,h ) = - i L 0 R ∂ t q i,h • ∂ t q i,h - i L 0 Γ ∂ 2 x,t q i,h • ∂ 2 x,t q i,h -a p d,h (∂ t U p,h , ∂ t U p,h ) - i r H Φ ′ H e i,h ∂ t e i,h 2 ≤ 0 (54)
In (54), which is the semi discrete equivalent of (29), the total energy of the system is given by

E tot (ξ h , q i,h , U p,h , U a,h ) = m H 2 ∂ t ξ h 2 + i k H Ψ H e i,h + i E s (q i,h ) + E p (U p,h ) + E a (U a,h ). ( 55
)
where E s,h (q i,h ) is the semi discrete energy of the i th string, E p,h (U p,h ) is the energy of the soundboard while E a,h (U a,h ) is the acoustic energy:

               E s,h (q i ) = 1 2 m s h (∂ t q i,h , ∂ t q i,h ) + L 0 H i (∂ x q i,h , q i,h ) dx (56a) E p,h (U p,h ) = 1 2 m p h (∂ t U p,h , ∂ t U p,h ) + 1 2 a p h (U p,h , U p,h ) (56b) E a,h (U a,h ) = 1 2 m a h (U a,h , U a,h ) (56c)
Theorem 1. The Cauchy problem associated with (42), or equivalently (53a), admits a unique global solution :

ξ h ∈ C 2 (R + ), q i,h ∈ C 2 (R + ; V s,∆x ), U p,h ∈ C 2 (R + ; V p,M,h ), U a,h ∈ C 2 (R + ; V a,h )
which moreover satisfies the stability result, for any T > 0: :

ξ h W 1,∞ (0,T ) ≤ C T, U a,h L ∞ (0,T,Ha) ≤ C T, q i,h W 1,∞ (0,T,Hs) + q i,h L ∞ (0,T,Vs) ≤ C T, U p,h W 1,∞ (0,T,Hp) + U p,h L ∞ (0,T,Vp) ≤ C T, (57)
where the constant C is independent of h and T .

Proof. It is a consequence of classical theorems on ordinary differential equations. One easily checks that the nonlinearity in (53a) (cf. in particular ( 5) and (8c)) is locally Lipschitz. This ensures the existence and uniqueness of a local solution in time. The fact that this solution is global is a consequence of the a priori estimates (57) that are themselves a consequence of the energy identity (54). Indeed, thanks to well-known properties of spectral finite elements, there exists a constant C 0 , independent of h, such that

m p h (∂ t q i,h , ∂ t q i,h ) ≥ C 0 ∂ t q i,h 2
Hs , m p h (∂ t U p,h , ∂ t U p,h ) ≥ C 0 ∂ t U p,h 2
Hp , a p h (U p,h , U p,h ) ≥ C 0 U p,h 2 Vp , m a h (U a,h , U a,h ) ≥ C 0 U p,h 2 Ha
while, on the other hand

L 0 H i (∂ x q i,h , q i,h ) dx ≥ C 0 q i,h 2 Vs
The details are left to the reader.

Practical considerations : choice of the space discretization parameters

Our goal is to make realistic computations on a Steinway model D which is roughly speaking 3 meters long and 2 meters large. One of course has to fix a choice for the discretization steps ∆x (for the string) and h (for the soundboard and for the air), the polynomial degrees r s , r p and r a and the number M of retained modes for the soundboard. As usual, such a choice is not only the result of a completely rigorous analysis. Many practical (the implementation issues are taken into account) and empirical considerations (trial and error procedures), that are difficult to explain in an article, have also influenced our choices. Nevertheless, is appeared to us to expose some of the reasons that have guided these choices, which is rarely done, also because this enters for a part in the justification of other choices that will be done for the time discretization.

Concerning the string, we wanted to realize an "optimal" compromise between accuracy and computational cost. For ensuring a precise enough approximation on the eigenfrequencies of the linearized model (see next section for a more rigorous definition) in the range of frequencies [0, 10] kHz, we have been led to choose :

∆x = L/200, r s = 4. ( 58
)
Remark 2. The specific frequency range [0, 10] kHz is of interest in this study for several reasons. First, it is included in the range of human audible frequencies [0, 20] kHz. Second, when complex sounds are concerned, that involve many frequencies, the human ear can only distinguish the frequency components which have the largest amplitude. To be more precise, a frequency content being 30 times less loud than the loudest one will not be heard.

Remark 3. It would be also possible to increase ∆x and increase accordingly r s , while keeping, for instance, the same number of degrees of freedom. This would probably improve the accuracy. However, this would also diminish the sparsity of the matrices A s h and R s h (and we do not speak of the nonlinear part here) and make "less local" the coupling at the bridge which would increase the computational cost. In this sense, (58) corresponds to a compromise.

The number of modes M was chosen in order that the set { 2πλ 1/2 m , m ≤ M } contains all the eigenfrequencies in the range [0, 10] kHz. This led us to take M = 2, 400.

(59) Concerning the 2D and 3D meshes, we wanted to use a mesh as close as possible to a regular mesh and the choice of the space step mesh h was first guided by the need to respect the geometry of the piano (in particular the rim) and the soundboard (in particular the ribs (represented in green in figure 3(a)) and the bridge (represented in magenta and yellow in figure 3(a))), which led us to choose for the soundboard This gives the mesh illustrated in figure 3. The 3D mesh in the air is obtained by extrusion of a 2D mesh delimiting the soundboard but also the rim and the acoustic region, as illustrated on figure 3(a), using a slightly larger stepsize in the z-direction (orthogonal to the soundboard, see figure 3(b)). Then, we choose the polynomial degrees r p and r m in such a way that

h = 2 cm. ( 60
)
• One guarantees an approximation of the eigenvalues λ m , m ≥ M with a satisfying relative error. For the particular case of the greater eigenvalue (for which the numerical error is maximal), a Q 8 approximation gives λ 2400 = 9, 909.7 Hz, while a Q 3 approximation gives λ 2400 = 10, 037.7 Hz (relative error of 1, 3 %) and a Q 4 approximation gives λ 2400 = 9921.4 Hz (relative error of 0, 12 %) Clearly, Q 3 approximation does not seem precise enough for our expectations while Q 4 does. This led us to

r p = 4. (61)
• As emphasized by dispersion analysis done in [START_REF] Fauqueux | Mixed finite elements with mass-lumping for the transient wave equation[END_REF], providing about 8 degrees of freedom by wavelength with 4 th polynomials leads to a very low numerical dispersion. This led us to

r a = 4. ( 62
)
4. Time discretisation

Generalities

The time discretization is the most tricky part of our numerical method. One reason is that it is not easy to use time discretization of arbitrary order. It is essentially based on a finite difference method with a constant time step ∆t. As usual, we shall use the superscript n as a time index, i.e. for any function of time u(t) :

u(t n ) ≃ u n , if t n = n∆t.
One of the reasons that requires a particular care to the time discretization is that it is not easy to use numerical schemes of arbitrary order, contrary to space discretization, so that we expect that the major numerical dispersion effects will be due to time stepping.

The discretization of the 3D acoustic part of the problem (see section 4.5) will rely on a standard explicit (this is mandatory for computational cost issues, due to the huge dimension N a h ) leap frog method on staggered grids (see section 4.2). The corresponding CFL condition will lead to our choice for the time step. Taking into account our choices for r a and h (see the previous section), this allows us to take as a time step ∆t = 10 -6 seconds.

(63)

On the contrary, as the number of the discrete string's unknowns is considerably smaller than the number of acoustic unknowns, we accept implicit schemes for the string's equations, which moreover appears more or less mandatory for stability reasons, because of the nonlinearity of the model. This is detailed in section 4.2.

Finally, for the time discretization of the plate equations, we shall use a specific procedure consisting in solving analytically the plate equations in each time interval [t n , t n+1 ], provided that the source terms coming from the strings and from the acoustic part of be model are frozen during this time interval. Here we exploit in an essential way the modal decomposition of the solution.This procedure was already used in the case of the guitar [START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF] : it can be done, as we shall see, in an unconditionally stable manner, which allows us to avoid the much too severe CFL limit that a standard explicit time stepping would impose.

Time discretization of the strings' equations

In this section, we are only interested in a single string which allows us to drop the index i, for the simplicity of notation, in q i , q i,h , Q i,h , T i 0 , H i , H i,2 , U i , A s i,h , etc ... This index will reappear in the next section devoted to a choir of chords.

Introduction

For the presentation, let us consider the time discretisation of the string's system alone, without any source term at the right hand side (the only data are thus initial conditions). Such a source will appear due to the coupling terms with the hammer and the soundboard.

M ∂ 2 t q + ∂ t R q -∂ x (Γ ∂ x q) -∂ x ∇ p H(∂ x q, q) + ∇ q H(∂ x q, q) = 0 (64)
that we rewrite as a perturbation of the corresponding linearized model, using :

M ∂ 2 t q + ∂ t R q -∂ x (Γ ∂ x q) + A s (∂ x ) q -∂ x U (∂ x q) = 0 (65)
where A s (∂ x ) is the second order (vector valued) differential operator defined by

A s (∂ x ) q := -∂ x (∇ p H 2 (∂ x q, q)) + ∇ q H 2 (∂ x q, q) (66)
In this section, we consider Dirichlet boundary conditions for (u i , v i ) and Neumann boundary conditions for ϕ i at both ends x = 0, L, which we can write in our notation:

Jq(x e , t) = 0, Id -J Γ ∂ 2 xt q + ∇ p H(∂ x q, q) (x e , t) = 0, x e = 0 or L (67)
The same finite element method as in section leads to the semi-discrete problem (in algebraic form)

M s h d 2 dt 2 Q h + R s h d dt Q h + A s h Q h + ∇U (Q h ) = 0 (68)
that we need to discretize in time.

Remark 4. In fact, the only difference with what we did for dealing with the strings equations in section 2 is that the space V s in (16) (and accordingly the finite element space V s,∆x in (33)), should be modified in order to account for the boundary conditions (67) :

V s = {q = (u, v, ϕ) ∈ H 1 (0, L) 2 / u(0) = v(0) = 0 ⇐⇒ Jq(0) = Jq(L) = 0}.
Let us first explain the time discretization of the linearized system:

M ∂ 2 t q + ∂ t R q -∂ x (Γ ∂ x q) -A s (∂ x ) q = 0 (69)
whose semi-discrete version is

M s h d 2 dt 2 Q h + R s h d dt Q h + A s h Q h = 0. ( 70
)
As already mentioned in section 4.1, one of the objective of the time discretization of (70) will be to ensure a weak numerical dispersion, namely the fact that the eigenfrequencies of the discrete problem should be very good approximations of eigenfrequencies of the continuous problem. These are defined in the dissipationless case without source term: R = Γ = 0, S = 0. These "continuous" eigenfrequencies are the values of f ∈ R for which there exists particular solution of the discrete problem of the form :

q(x, t) = q 0 (x) e 2iπf t with q 0 : [0, L] → C N , q 0 = 0, while the "discrete" eigenfrequencies are defined as the values of f ∈ R for which there exists particular solution of the fully discrete problem (to be constructed) with unknowns Q n h ∈ R N s h (that represents the vector of degrees of freedom on the approximate solution at time T n = n ∆t) of the form :

Q n h = Q h,0 e 2iπf t with Q h,0 = 0.
One sees in particular that f is a continuous eigenfrequency if and only if 4πf 2 is an eigenvalue of the operator A s (∂ x ) with boundary conditions (67). This is a selfadjoint eigenvalue problem whose solutions have been described in detail in the first part of the present article [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF]. The eigenfrequencies can be split into three parts :

f n , n ≥ 1 ∩ f ℓ n , n ≥ 1 ∩ f S n , n ≥ 1 (71) where
• f n , n ≥ 1 are the eigenfrequencies associated with the flexural modes,

• f ℓ n , n ≥ 1 are the eigenfrequencies associated with the longitudinal modes , • f S n , n ≥ 1 are the eigenfrequencies associated with the shear modes.

In the case of the piano, we have to privilege the approximation of the flexural modes for the following reasons:

• The flexural modes are the ones that are priorly solicited by the hammer ,

• These are the modes that define the note played by the piano,

• Most eigenfrequencies located in the range of audible eigenfrequencies are associated to flexural modes.

Roughly speaking, a good approximation of flexural modes should be achieved if we perform an accurate approximation of the part of the model that corresponds to flexural modes. Of course, in the case of a coupled system, this is not a very well defined part. However, it is known that, in the small radius limit, the linear model can be seen as an approximation of the D'Alembert's model (see section 3.2.1 of [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF]). This corresponds to the following decomposition of the quadratic function (the superscript D here refers to D'Alembert while the superscript p refers to "perturbation")

H 2 (p, q) = H D 2 (p, q) + H p 2 (p, q), H D 2 (p, q) = 1 2 A D p • p, A D =   T 0 0 0 0 0 0 0 0 0   . ( 72
)
The reader will easily check that (see also [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF], formula (60))

H p 2 (p, q) = EA 2 |p 2 | 2 + |p 3 | 2 + AGκ 2 |q 3 -p 1 | 2 + (EA -T 0 ) 1 2 |p 1 | 2 + (1 + p 2 ) -p 2 1 + (1 + p 2 ) 2 ≥ 0. ( 73
)
This corresponds to an analogous decomposition of the operator A s (∂ x ) and the corresponding matrix A s h :

A s (∂ x ) = A s,D (∂ x ) + A s,p (∂ x ), A s h = A s,D h + A s,p h (74)
where, according to (73), the selfadjoint operator A s,p (∂ x ) and the corresponding symmetric matrix A s,p h are both positive. In other words, we see (69) as a perturbation of the simpler model

M ∂ 2 t q + ∂ t R q -∂ x (Γ ∂ x q) + A s,D (∂ x ) q = 0 (75)
with semi-discrete version

d 2 dt 2 M s h Q h + d dt R s h Q h + A s,D h Q h = 0 (76)
As explained in [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF], due to the smallness of the diameter of the strings, the eigenfrequencies of (75) are good approximations of the the eigenfrequencies of the flexural modes, at least in the range of audible frequencies (the difference representing precisely the inharmonicity effects).

Construction of the scheme.

According to the previous paragraph, the idea is thus to start from an accurate time discretisation of (76). It is well known that the θ-scheme exhibits a super convergence for θ = 1/12, leading to fourth order accuracy (see [START_REF] Chabassier | Introduction and study of fourth order theta schemes for linear wave equations[END_REF]). That is why we choose to start with the following discretization:

M s h Q n+1 h -2 Q n h + Q n-1 h ∆t 2 + R s h Q n+1 h -Q n-1 h 2∆t + A s,D h {Q h } n 1 12 = 0 (77)
where by definition:

{Q h } n θ := θ Q n+1 h + (1 -2 θ) Q n h + θ Q n-1 h . ( 78
)
The negative counterpart of the fourth order accuracy of ( 77) is that it is submitted to a CFL stability condition

∆t ρ sp (M s h ) -1 A s,D h ≤ 3 2 . ( 79
)
where ρ sp (•) holds for the spectral radius of a matrix. This gives a condition of the form

c f ∆t ∆x ≤ C(r s ) (80)
where

c f = (T 0 /ρ) 1 2
is the velocity of flexural waves and C(r s ) a constant that depends only on the polynomial degree r s . However, when ∆x and r s are given by (58), this constraint is compatible with the time step (63) imposed by the discretization of the 3D acoustic equation.

It is a priori tempting to apply the same 1/12-scheme to the full linearized system (70) but we had to abandon this idea for stability reasons. Indeed, due again to the high velocity contrast between longitudinal or shear waves versus flexural waves, the matrix (M s h ) -1 A s h appears to have a much higher spectral radius than (M s h ) -1 A s,D h : the stability condition is roughly speaking of the form (80) where c f is replaced by the maximal velocity between shear and longitudinal waves, which is much higher than c f . As a consequence the stability condition would lead to a much smaller time step. That is why we have chosen to approximate the "perturbation part" of the equation, in a less accurate (second order) but unconditionally stable way. For this, choosing θ ≥ 1/4, we consider the following scheme, according to (72, 74) (a particular case of the (θ, θ)-schemes studied in [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF]):

M s h Q n+1 h -2Q n h + Q n-1 h ∆t 2 + R s h Q n+1 h -Q n-1 h 2∆t + A s,D h {Q h } n 1 12 + A s,p h {Q h } n θ = 0 (81)
Finally, to construct our numerical scheme for (68), it remains to include the nonlinear term ∇U (Q h ). One major difficulty was to preserve the stability of the scheme and, if possible, discretize in time the non linear term in such a way that the stability condition (79) is not affected. This is obtained by constructing an "energy preserving scheme". For this, we rely on the specific developments we made in [START_REF] Chabassier | Energy preserving schemes for nonlinear hamiltonian systems of wave equations. application to the vibrating piano string[END_REF] for constructing conservative discretizations of hamiltonian nonlinear systems of wave equation. The scheme proposed in [START_REF] Chabassier | Energy preserving schemes for nonlinear hamiltonian systems of wave equations. application to the vibrating piano string[END_REF] can be seen alternatively as an extension to systems of the Strauss-Vasquez scheme for scalar equations (see [START_REF] Strauss | Numerical solution of a nonlinear klein-gordon equation[END_REF]), or an extension to nonlinear problems of the well-known θ-schemes (or conservative Newmark schemes, see [DLB + 00]) for linear problems. Note that other energy preserving finite difference schemes for piano strings have been proposed by [START_REF] Bilbao | Conservative numerical methods for nonlinear strings[END_REF] in the case of polynomial nonlinearities. Our scheme is based on the notion of conservative approximate gradient that we recall now :

∇U (Q h )(t n ) is approximated by ∇U (Q h )(t n ) ∼ ∇h U (Q n+ 1 2 h , Q n-1 2 h ) with Q n+ 1 2 h = Q n+1 h + Q n h 2 , (82)
and where ∇h U : R N h × R N h → R N h is the so-called approximate gradient (see also [START_REF] Quispel | Discrete gradient methods for solving odes numerically while preserving a first integral[END_REF]) that is built in order to ensure the consistency of the approximation (82) as well as the conservativity property:

∇h U (Q + h , Q - h ) • Q + h -Q - h = U h (Q + h ) -U h (Q - h ), (83)
where the discrete energy density U h (Q h ) is defined by:

U h (Q h ) • Q h := L 0 U ∂ x q h ∂ x q h ( 84
)
and q h (resp. q h ) is the element of V s,∆x associated with the vector of degrees of freedom Q h (resp. Q h ).

Remark 5. The important consequence of (83) is the identity:

∇h U Q n+ 1 2 h , Q n-1 2 h • Q n+ 1 2 h -Q n-1 2 h 2∆t = 1 ∆t U h Q n+ 1 2 h -U h Q n-1 2 h
which appears as a discrete equivalent of the chain rule

d dt L 0 U (∂ x q) = L 0 ∇U (∂ x q) • ∂ 2 xt q
that is used for the proof of the energy identity for the continuous problem.

The definition of ∇h U requires preliminary notation.

Notation : Let Φ : R N → R such that Φ(u 1 , . . . , u N ) ∈ R. We set Φ ℓ (u ℓ ; u j =ℓ ) := Φ(u 1 , . . . , u N ) and

∀ (x, x ′ ) ∈ R 2 , ∀ y ∈ R N -1 , δ ℓ Φ(x, x ′ ; y) =      Φ ℓ (x ; y) -Φ ℓ (x ′ ; y) x -x ′ if x = x ′ ∂ ℓ Φ ℓ (x ; y) if x = x ′ (85)
where ∂ ℓ Φ ℓ (x ; y) := ∂Φ ∂u ℓ (x; y) with an obvious abuse of notation.

In the sequel, N = 3 and for simplicity we denote q ℓ , 1 ≤ ℓ ≤ 3 the ℓ th component of q (note that, according to the notation of section S, q 1 ≡ u, q 2 ≡ v, q 3 ≡ ϕ ) and, for avoiding two many indices, we use the same notation for the components of q h . Finally, we denote P the set of permutations of {1, 2, 3} with #P = 6.

Then, for any

(Q + h , Q - h ) ∈ R N h × R N h , we define ∇h U (Q + h , Q - h ) ∈ R N h such that, for any Q h ∈ R N h ∇h U (Q + h , Q - h ) • Q h := 1 6 L 0 π∈P 3 ℓ=1 δ ℓ U (∂ x q + ℓ , ∂ x q - ℓ ; ∂ x q sπ(ℓ,j) j =ℓ ) • ∂ x q ℓ , s π (ℓ, j) = sgn π(j) -π(ℓ) (86)
where q + h (resp. q - h , q h ) is the element of V associated with the vector

Q + h (resp. Q - h , Q h ).
The reader can check that, as it is proven in a more general case in [START_REF] Chabassier | Energy preserving schemes for nonlinear hamiltonian systems of wave equations. application to the vibrating piano string[END_REF], that the function ∇h U defined by (86) does satisfy the conservativity property (83).

Finally, according to (82), the numerical scheme that we propose for the time discretization of (68) is

M s h Q n+1 h -2Q n h + Q n-1 h ∆t 2 + R s h Q n+1 h -Q n-1 h 2∆t + A s,D h {Q h } n 1 12 + A s,p h {Q h } n θ + ∇h U (Q n+ 1 2 h , Q n-1 2 h ) = 0 (87)
It is not difficult (see again [START_REF] Chabassier | Energy preserving schemes for nonlinear hamiltonian systems of wave equations. application to the vibrating piano string[END_REF]) to check theoretically and numerically that this scheme is second order accurate in time.

Stability analysis

It is possible to prove a stability result via a discrete energy identity. This is the object of the following Proposition 1. Any solution to (81) satisfies the following energy identity:

E n+ 1 2 s,h -E n-1 2 s,h ∆t = -R s h Q n+1 h -Q n-1 h 2∆t • Q n+1 h -Q n-1 h 2∆t ≤ 0 ( 88
)
where the discrete string's energy is defined by:

E n+ 1 2 s,h = 1 2 M s,θ h Q n+1 h -Q n h ∆t • Q n+1 h -Q n h ∆t + 1 2 A s h Q n+ 1 2 h • Q n+ 1 2 h + U h Q n+ 1 2 h ( 89
)
where the modified mass matrix M s,θ h is given by:

M s,θ h := M s h + θ - 1 4 ∆t 2 A s,p h - ∆t 2 6 A s,D h (90)
Proof. We use the well-known identity

{Q h } n θ = {Q h } n 1 4 + θ - 1 4 ∆t 2 Q n+1 h -2 Q n h + Q n-1 h ∆t 2
, to rewrite (87) as

M s,θ h Q n+1 h -2Q n h + Q n-1 h ∆t 2 + R s h Q n+1 h -Q n-1 h 2∆t + A s h {Q h } n 1 4 + ∇h U (Q n+ 1 2 h , Q n-1 2 h ) = 0. (91)
We take the scalar product of (77) with

Q n+1 h -Q n-1 h 2∆t ≡ Q n+ 1 2 h -Q n-1 2 h ∆t ,
use the the conservativity property (83) and the well known identities:

M s,θ h Q n+1 h -2Q n h + Q n-1 h ∆t 2 • Q n+1 h -Q n-1 h 2∆t = 1 2∆t M s,θ h Q n+1 h -Q n h ∆t • Q n+1 h -Q n h ∆t -M s,θ h Q n h -Q n-1 h ∆t • Q n h -Q n-1 h ∆t , A s h {Q h } n 1 4 • Q n+1 h -Q n-1 h 2∆t = 1 2∆t A s h Q n+ 1 2 h • Q n+ 1 2 h -A s h Q n-1 2 h • Q n-1 2 h .
We then easily obtain the announced result.

Corollary 1. The numerical scheme (77) is stable as soon as θ ≥ 1 4 and the stability condition (79) is satisfied.

Proof. We shall not enter all the details which are more or less standard manipulation (see also remark 6). The stability result appears as a consequence of the positivity of the energy (89). According to the definition of the matrix A s h (43) and of the function U h (84), we observe that

1 2 A s h Q n+ 1 2 h • Q n+ 1 2 h + U h Q n+ 1 2 h = L 0 H 2 (∂ x q n+ 1 2 h , q n+ 1 2 h ) + U (∂ x q n+ 1 2 h ) dx = L 0 H(∂ x q n+ 1 2 h , q n+ 1 2 h ) dx ≥ 0,
because of the positivity of H and the positivity of the weights in the Gauss Lobatto formulas. Therefore, the positivity of ( 89) is a consequence of the positivity of the modified mass matrix the modified mass matrix M s,θ h . This matrix is in fact the sum of two positive matrices since

• θ - 1 4 ∆t 2 A s,p
h is positive because of θ ≥ 1/4 and (73),

• M s h - ∆t 2 6 A s,D h is positive thanks to (79).
This concludes the proof.

Remark 6. In corollary 1, we speak of stability in a rather vague way. In fact, proceeding as in [START_REF] Joly | Variational methods for time-dependent wave propagation problems[END_REF] for instance, see also [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF], it is possible, under suitable assumptions on the source term, to establish uniform estimates in discrete L ∞ 0, T ; H 1 (0, L) 3 norm :

sup t n ≤T q n h H 1 (0,L) ≤ C T.
Remark 7. Since the problem is nonlinear, the stability result of corollary 1 does not imply the existence and uniqueness of the discrete solution. At each time step, to find Q n+1 h one has to solve a nonlinear equation of the form

M s h + 2∆t R s h + ∆t 2 A s h Q n+1 h + ∆t 2 F n h Q n+1 h = (rhs) n h (92)
where (rhs) n h is a known right hand side (from the two previous instants) and the nonlinear function F n h : R N → R N is defined by

F n h Q h := ∇h U Q h + Q n-1 h 2 , Q n-1 2 h ( 93
)
We have not investigated in detail the (unique) solvability of equation (92), which is not necessarily a trivial issue (this is a weak point of our analysis). We conjecture that, exploiting the invertibility of M s h + 2∆t R s h + ∆t 2 A s h under the condition (79), in other words considering the nonlinearity as a perturbation, we should be able to prove by a fixed point method the existence and uniqueness of the solution at least for ∆t/∆x small enough.

Moreover, in practice, we never met any serious difficulty in the numerical resolution of (92).

Time discretization of the strings' choir / hammer system

This section is devoted to the time discretisation of the coupling between the strings' choir and a nonlinear hammer. This corresponds to the following subsystem of (53a), in which we have removed at the right hand side of the string equation the term F p i,h (t) L p,h + F r i,h (t) L r,h that is due to the coupling with the soundboard:

                           Find ξ h , (Q i,h ) 1≤i≤Ns : R + → R × (V s,∆x ) Ns (94a)      m H d 2 ξ h dt 2 (t) = - i F i,h (t), F i,h (t) = k H Φ H e i,h (t) + r H d dt Φ H e i,h (t) , e i,h (t) = ξ + ξ i,h (t) -ξ h (t) + , ξ i,h (t) = L H,h • Q i,h (94b)
M s h d 2 dt 2 Q i,h + R s h d dt Q i,h + A s i,h Q i,h + ∇U i (Q i,h ) = F i,h (t) L H,h , ∀ 1 ≤ i ≤ N s (94c)
The numerical scheme we propose is based on equation (87) for the string equation and is inspired by [START_REF] Rhaouti | Time-domain modeling and numerical simulation of a kettledrum[END_REF] for the hammer's equation. The main point is the approximation of the source term F i,h (t). Introducing Φ + H (ν) = Φ H (ν + ) and a primitive Ψ + H of Φ + H , we remark that

F i,h (t) = k H Ψ + H ′ d i,h (t) + r H d dt Φ + H d i,h (t)
where

d k i,h = ξ + ξ k i,h -ξ k h such that e k i,h = (d k i,h ) + . Let us introduce the application δΨ + H (a, b) = Ψ + H (a) -Ψ + H (b) a -b , if a = b, δΨ + H (a, a) = Ψ ′ H (a) ≡ Φ H (a) (95)
we use the second order approximation :

F i,h (t n ) ∼ k H δΨ + H d n+1 i,h , d n-1 i,h + r H Φ + H (d n+1 i,h ) -Φ + H (d n-1 i,h ) 2∆t which leads to the sheme                            m H ξ n+1 h -2ξ n h + ξ n-1 h ∆t 2 = - i F n i,h , F n i,h = -k H δΨ + H d n+1 i,h , d n-1 i,h + r H Φ + H (d n+1 i,h ) -Φ + H (d n-1 i,h ) 2∆t , d k i,h = ξ + ξ k i,h -ξ k h , ξ k i,h = L H,h • Q k i,h , k = n -1, n + 1 M s h Q n+1 i,h -2Q n i,h + Q n-1 i,h ∆t 2 + R s h Q n+1 i,h -Q n-1 i,h 2∆t + A s,D i,h {Q i,h } n θ + A s,p i,h {Q i,h } n 1 12 + ∇U (Q n+ 1 2 i,h , Q n-1 2 i,h ) = F n i,h L H,h (96)
where we recall that Q

n± 1 2 i,h := Q n i,h + Q n±1 i,h 2
. This system is completed with initial data

ξ 0 h = -ξ, ξ 1 h = -ξ + v H ∆t, Q 0 i,h = Q 1 i,h = 0, ∀ 1 ≤ i ≤ N s . ( 97
)
The stability of (96) then follows again from a discrete energy estimate :

Proposition 2. Any solution to the scheme (96) satisfies the following energy identity:

E n+ 1 2 H,s,h -E n-1 2 H,s,h ∆t = - i k H Φ + H (d n+1 i,h ) -Φ + H (d n-1 i,h ) 2∆t d n+1 i,h -d n-1 i,h 2∆t - i R s h Q n+1 i,h -Q n-1 i,h 2∆t • Q n+1 i,h -Q n-1 i,h 2∆t (98)
where the discrete energy E n+ 1 2 H,s,h is defined by

E n+ 1 2 H,s,h = i E n+ 1 2 i,s,h + E n+ 1 2 H,h , E n+ 1 2 H,h := m H 2 ξ n+1 h -ξ n h ∆t 2 + k H Ψ H (e n+1 i,h ) + Ψ H (e n i,h ) 2 (99)
where E n+ 1 2 i,s,h is the discrete energy associated with i th string, defined by (89) (having reincorporated the index i).

Proof. It is quite standard. We multiply the first equation in ( 96) by

ξ n+1 h -ξ n-1 h /2∆t to obtain m H 2∆t ξ n+1 h -ξ n h ∆t 2 - ξ n h -ξ n-1 h ∆t 2 = - i F n i,h ξ n+1 h -ξ n-1 h 2∆t .
On the other hand, proceeding as in the proof of proposition 1, we get

E n+ 1 2 i,s,h -E n-1 2 i,s,h ∆t = F n i,h L H,h • Q n+1 h -Q n-1 h 2∆t -R s h Q n+1 h -Q n-1 h 2∆t • Q n+1 h -Q n-1 h 2∆t .
After summation of these equalities, we get

1 ∆t m H 2 ξ n+1 h -ξ n h ∆t 2 + i E n+ 1 2 i,s,h - m H 2 ξ n h -ξ n-1 h ∆t 2 + i E n-1 2 i,s,h + i F n i,h d n+1 i,h -d n-1 i,h 2∆t = 0.
Finally, using the expression of F n i,h and the definition of δΨ + H , we have

F n i,h d n+1 i,h -d n-1 i,h 2∆t = k H Ψ + H d n+1 i,h -Ψ + H d n-1 i,h 2 + r H Φ + H d n+1 i,h -Φ + H d n-1 i,h 2∆t
d n+1 i,h -d n-1 i,h 2∆t (100)
Substituting (100) in the previous equality leads to (98), using (99) and the link between Φ + H and Φ H .

Corollary 2. Under the same assumptions as in corollary 1, the scheme (96) is stable.

Proof. The matrix R s h being positive and the function Φ + H increasing, (98) means that the energy E

n+ 1 2 H,s,h is a decreasing sequence. It is thus bounded by its initial value E n+ 1 2 H,s,h = E 1 2 H,s,h ≡ 1 2 |v H | 2
which provides a uniform bound with respect to ∆t and h. The assumptions of corollary 1, in particular the stability condition (79), imply the positivity of the energy and thus (see corollary 1 and remark 6) the stability of the scheme.

Time discretization of the strings / soundboard system

This section is devoted to the time discretisation of the coupling between the strings' choir and the soundboard. This corresponds to the following subsystem of (53a), in which we have removed the term F i,h (t) L H,h at the right hand side of the string equation that is due to the coupling with the hammer and the term t C ap,h U a,h at the right hand side of the soundboard equation that is due to the coupling with the acoustics:

                             Find (Q i,h ) 1≤i≤Ns , U p,h , F p i,h , F r i,h : R + → R N s h × R N p h × R × R (101a) M s h d 2 dt 2 Q i,h + R s h d dt Q i,h + A s i,h Q i,h + ∇U i (Q i,h ) = F p i,h (t) L p,h + F r i,h (t) L p,h ∀ 1 ≤ i ≤ N s (101b) L p,h • Q i,h = L s,h • U p,h , L r,h • Q i,h = 0 (101c) d 2 dt 2 U p,h + Λ p d,h d dt U p,h + Λ p h U p,h = - i F p i,h (t) L s,h (101d)
Of course, we shall rely on the scheme (87) for the time discretization of (101b) : only the way we discretize F p i,h (t) and F p i,h (t) needs to be clarified. On the other hand, we could also use a finite difference scheme for the soundboard equation 101d. We have abandoned this idea to avoid the introduction of an additional source of error due to numerical dispersion : we opted for a quasi-exact resolution that moreover exploits our choice of a modal basis and the diagonal nature of the matrices Λ p h and Λ p d,h . This strategy was already adopted for the simulation of the guitar [START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF] and it permits to avoid any time step limitation related to stability condition (see [START_REF] Cohen | Mixed higher order spectral finite elements for reissner-mindlin equations[END_REF]) and to avoid alternative local time stepping procedures (see [START_REF] Diaz | Energy conserving explicit local time-stepping for second-order wave equations[END_REF][START_REF] Fouquet | Conservative space-time mesh refinement methods for the fdtd solution of maxwell's equations[END_REF]).

Time discretization of the soundboard equation

We assume that the right hand side F p i,h (t) is approached by a piecewise constant function on intervals [t n-1 2 , t n+ 1 2 ]:

F p i,h (t) ≃ F p,n i,h , for t ∈ [t n-1 2 , t n+ 1 2 ] ( 102
)
The simple idea is then to compute the exact solution associated with this discretized right hand side. As a result, the discrete unknowns for U p,h will be

U n+ 1 2 p,h ≃ U p,h (t n+ 1 2 ), Un+ 1 2 p,h ≃ dU p,h dt (t n+ 1 2 ) (103)
and the passage from t n-1 2 to t n+ 1 2 will be given by:

U n+ 1 2 p,h = U n p,h (t n+ 1 2 ), Un+ 1 2 p,h = d U n p,h dt (t n+ 1 2 ) ( 104
)
where

U n p,h (t) : t n-1 2 , t n+ 1 2 → R N a,h is the solution of d 2 dt 2 U n p,h (t) + Λ p d,h d dt U n p,h (t) + Λ p h U n p,h (t) = - i F p,n i,h . L s,h . (105)
Solving explicitly (105) leads to the formulas :

       U n+ 1 2 p,h = S 0 ∆t U n-1 2 p,h + S 1 ∆t Un-1 2 p,h - i F p,n i,h R ∆t L s,h Un+ 1 2 p,h = Ṡ0 ∆t U n-1 2 p,h + Ṡ1 ∆t Un-1 2 p,h - i F p,n i,h Ṙ∆t L s,h (106)
where (S 0 ∆t , S 1 ∆t , R ∆t ) and ( Ṡ0 ∆t , Ṡ1 ∆t , Ṙ∆t ) are N p h ×N p h diagonal matrices whose respective entries (s m,0 ∆t , s m,1 ∆t , r m ∆t ) and ( ṡm,0 ∆t , ṡm,1 ∆t , ṙm ∆t ) are given by the explicit formulas

               s m,0 ∆t = e -α h,m ∆t λh,m cos λm,h ∆t + α h,m λm,h sin( λm,h ∆t) s m,1 ∆t = e -α h,m ∆t λh,m sin( λm,h ∆t) r m ∆t = 1 λ h,m λh,m λ m,h 1 -e -α h,m ∆t cos λm,h ∆t -α h,m e -α h,m ∆t sin( λm,h ∆t) (107)                ṡm,0 ∆t = e -α h,m ∆t λh,m cos λm,h ∆t + α h,m λm,h sin( λm,h ∆t) ṡm,1 ∆t = e -α h,m ∆t λh,m sin( λm,h ∆t) ṙm ∆t = 1 λ h,m λh,m λ m,h 1 -e -α h,m ∆t cos λm,h ∆t -α h,m e -α h,m ∆t sin( λm,h ∆t) (108)
where we have defined

α m,h = 1 2 f ve (λ h,m ) ∈ R, λm,h = λ h,m -α 2 m,h 1 2 ∈ R ∪ i R
and by convention, for any real x, x

1 2 = √ x if x ≥ 0 and x 1 2 = i |x| if x ≤ 0.

Time discretisation of the coupled problem

According to the previous section, assuming that F r i,h (t) is approached by a piecewise constant function as F p i,h (t) (see (102)), the time discretization of equation (101b) is straightforward :

M s h Q n+1 i,h -2Q n i,h + Q n-1 i,h ∆t 2 + R s h Q n+1 i,h -Q n-1 i,h 2∆t + A s,D i,h {Q i,h } n θ + A s,p i,h {Q i,h } n 1 12 + ∇U (Q n+ 1 2 i,h , Q n-1 2 i,h ) = F p,n i,h L p,h + F r,n i,h L r,h (109)
To complete (109) and (106), we need to discretize the coupling equation (101c), which is not that obvious.

Once again, we shall do it in an energy preserving (or energy decaying) way to ensure stability.

Defining the discrete soundboard energy as

E n+ 1 2 p,h = 1 2 Un+ 1 2 p,h • Un+ 1 2 p,h + 1 2 Λ p h U n+ 1 2 p,h • U n+ 1 2 p,h (110)
we observe that it satisfies the identity:

E n+ 1 2 p,h -E n-1 2 p,h ∆t + t n+ 1 2 t n-1 2 Λ p d,h d dt U p,h (s) • d dt U p,h (s) ds = - i F p,n i,h L s,h • U n+ 1 2 p,h -U n-1 2 p,h ∆t (111)
To derive this identity, it suffices to take the inner product of (105) with d dt U p,h , and integrate the resulting equality in time between t n -∆t/2 and t n + ∆t/2.

On the other hand, taking the inner product of (109) by

Q n+1 i,h -Q n-1 i,h 2∆t
and proceeding as for proposition 1, we get

E n+ 1 2 i,s,h -E n-1 2 i,s,h ∆t + R s h Q n+1 i,h -Q n-1 i,h 2∆t • Q n+1 i,h -Q n-1 i,h 2∆t = F p,n i,h L p,h + F r,n i,h L r,h • Q n+1 i,h -Q n-1 i,h 2∆t .
Therefore, defining the discrete string-plate energy as

E n+ 1 2 s,p,h = E n+ 1 2 p,h + i E n+ 1 2 i,s,h (112)
and imposing, from the second coupling equation of (101c),

L r,h • Q n i,h = 0 (113)
we get from ( 112) and (111) the energy identity

E n+ 1 2 s,p,h -E n-1 2 s,p,h ∆t + i R s h Q n+1 i,h -Q n-1 i,h 2∆t • Q n+1 i,h -Q n-1 i,h 2∆t + t n+ 1 2 t n-1 2 Λ p d,h ∂ t U p,h (s) • ∂ t U p,h (s)
= i F p,n i,h L p,h • Q n+1 i,h -Q n-1 i,h 2∆t -F p,n i,h L s,h • U n+ 1 2 p,h -U n-1 2 p,h
∆t .

To ensure an appropriate energy transfer between the strings and the soundboard, it is natural to impose

L p,h • Q n+1 i,h -Q n-1 i,h 2∆t = L s,h • U n+ 1 2 p,h -U n-1 2 p,h ∆t (114)
which appears as a second order accurate time discretization of the first coupling equation of (101c) after time differentiation, i.e. in the form

L p,h • dQ i,h dt = L s,h • dU p,h dt .
Finally, the fully discrete scheme we propose for the approximation of is made of equations (105,109,113,114).

From what precedes, we deduce that for any solution of this scheme we have the identity

E n+ 1 2 s,p,h -E n-1 2 s,p,h ∆t + i R s h Q n+1 i,h -Q n-1 i,h 2∆t • Q n+1 i,h -Q n-1 i,h 2∆t + t n+ 1 2 t n-1 2 Λ p d,h d dt U p,h (s) • d dt U p,h (s) = 0 (115)
which corresponds to a discrete energy decay result, thanks to the positivity of the matrices R s h and Λ p d,h . Proceeding as in the previous section, it is then easy to show the following stability result: Proposition 3. Under the same assumptions as in corollary 1, the scheme (105, 109, 113, 114) is stable.

Discretization of the soundboard-acoustics system

In this section, we are interested in the time discretization of the fluid-structure interaction problem describing the coupling between acoustics and soundboard vibrations. This corresponds to the following subsystem made of the last two equations of (53a), namely (53e) and (53f), in which we have simply removed from the right hand side of equation (53e) the source term due to the coupling with the strings, i. e. assuming that

F p i,h (t) ≡ 0 :        d 2 dt 2 U p,h + Λ p d,h d dt U p,h + Λ p h U p,h = C t ap,h U a,h M a h d dt U a,h + B a h U a,h = -C ap,h d dt U p,h (116)

Time discretisation of the acoustic propagation

We first consider the time discretization of (here S h (t) denotes any discretized source term)

M a h d dt U a,h + B a h U a,h = S h (t) (117)
Given the size of the computational domain for the acoustic propagation (more than 90 000 000 degrees of freedom), an explicit strategy is needed. We have chosen to use a classical interleaved version of the leap frog scheme for the first order system (118), in which the two components of the discrete acoustic unknown vector U a,h , namely the discrete pressure P h and the discrete velocity V h are evaluated on two staggered grids in time : P h at times t n+ 1 2 and V h at times t n , a choice which will justified when studying the stability of the vibroacoustic coupling. For this, assuming for simplicity that the source term appears only in the sound pressure equation, and using remark 1, more precisely (51), we rewrite (117) as

     M p h d dt P h + C h V h = S h (t) M V h d dt V h -C t h P h = 0 (118)
which leads to the scheme

           M p h P n+ 1 2 h -P n-1 2 h ∆t + C h V n h = S n h (119a) M V h V n+1 h -V n h ∆t -C t h P n+ 1 2 h = 0 (119b)
The stability analysis of this scheme is very classical. We recall it for completeness.

Proposition 4. Any solution to (119) satisfies the energy identity:

E n+ 1 2 a,h -E n-1 2 a,h ∆t = S n h • P n+ 1 2 h + P n-1 2 h 2 (120)
where the acoustical energy is defined as

E n+ 1 2 a,h = 1 2 M p h P n+ 1 2 h • P n+ 1 2 h + 1 2 M V h V n+1 h • V n h ( 121
)
As a consequence, the scheme is L 2 -stable as soon as ∆t is chosen small enough so that

the matrix M p h - ∆t 2 4 C t h (M V h ) -1 C h is positive. (122)
Proof. Obtaining (121) is left to the reader. Next, to get (122), it suffices to observe

E n+ 1 2 a,h = 1 2 M p h P n+ 1 2 h • P n+ 1 2 h + 1 2 M V h V n+1 h + V n h 2 • V n+1 h + V n h 2 - ∆t 2 8 C t h (M V h ) -1 C h P n+ 1 2 h • P n+ 1 2 h
Remark 8. In practice, condition (122) implies an upper bound on the time step ∆t as for condition (79).

Time discretization of the coupled system

In this section, we go back to the time discretization of the system (116). According to the previous paragraph, this system can be rewritten using the decomposition of the acoustic unknown U a,h as U a,h ≡ (P h , V h ) t . To do so, we use (52) to rewrite (116) as:

                 d 2 dt 2 U p,h + Λ p d,h d dt U p,h + Λ p h U p,h = C t pp,h P h      M p h d dt P h + C h V h = -C pp,h d dt U p,h M V h d dt V h -C t h P h = 0 (123)
Our scheme will be inspired from what has been done in section 4.4.1 for the plate unknowns, here for the first equation of (123), and from section 4.5.1 for the acoustic unknowns. In particular, the discrete unknowns are:

U n+ 1 2 p,h ≃ U p,h t n+ 1 2 , Un+ 1 2 p,h ≃ dU p,h dt t n+ 1 2 P n+ 1 2 h ≃ P h t n+ 1 2 , V n h ≃ V h t n (124)
The specificity of our coupling scheme holds in the way we discretize / approximate the right hand sides in (123), which we do in a particular centered manner. More precisely • for the first equation of (123), during the time interval [ t n-1 2 , t n+ 1 2 ], we solve (123) exactly, assuming that C t pp,h P h is approximated by a constant in time:

C t pp,h P h ≃ C t pp,h P n+ 1 2 h + P n-1 2 h 2 in [ t n-1 2 , t n+ 1 2 ]
Let us recall that, in practice, this can be done explicitly with formulas analogous to (106) where

- i F p,n i,h R ∆t L s,h is replaced by C t pp,h P n+ 1 2 h + P n-1 2 h 2
• for the first equation of (123), we use the numerical scheme (119) with

S n h := -C t pp,h U n+ 1 2 h -U n-1 2 h ∆t as a approximation of -C t pp,h d dt U h (t n ).
This results into the following scheme

                                               U n+ 1 2 p,h = U n p,h (t n+ 1 2 ), Un+ 1 2 p,h = d U n p,h dt (t n+ 1 2 ) where (125a)          d 2 dt 2 U n p,h (t) + Λ p d,h d dt U n p,h (t) + Λ p h U n p,h (t) = C t pp,h P n+ 1 2 h + P n-1 2 h 2 , in t n-1 2 , t n+ 1 2 U n p,h (t n-1 2 ) = U n-1 2 p,h , d dt U n p,h (t n-1 2 ) = Un-1 2 p,h (125b)
         M p h P n+ 1 2 h -P n-1 2 h ∆t + C h V n h = -C pp,h U n+ 1 2 p,h -U n-1 2 p,h ∆t M V h V n+1 h -V n h ∆t -C t h P n+ 1 2 h = 0 (125c)
Proposition 5. Any solution to (125) satisfies the following energy identity:

E n+ 1 2 p,a,h -E n-1 2 p,a,h ∆t = - t n+ 1 2 t n-1 2 Λ p d,h d dt Ũn p,h (s) • d dt Ũn p,h (s) ds (126)
where the energy of the vibro-acoustic system is

E n+ 1 2 p,a,h = E n+ 1 2 p,h + E n+ 1 2 a,h (127)
where E n+ 1 2 p,h is defined by (110) and E n+ 1 2 a,h by (121). As a consequence, under assumption (122), the numerical scheme (125) is stable.

Proof. Proceeding as in section 4.4.1, we get the following energy identity for the soundboard unknowns:

E n+ 1 2 p,h -E n-1 2 p,h ∆t = C t pp,h P n+ 1 2 h + P n-1 2 h 2 • U n+ 1 2 p,h -U n-1 2 p,h ∆t - t n+ 1 2 t n-1 2 Λ p d,h d dt Ũn p,h (s) • d dt Ũn p,h (s) ds. (128)
On the other hand, using proposition 4, we get the following energy identity for the acoustic variables:

E n+ 1 2 a,h -E n-1 2 a,h ∆t = -C pp,h U n+ 1 2 p,h -U n-1 2 p,h ∆t • P n+ 1 2 h + P n-1 2 h 2 (129)
and adding (128) and (129), the coupling terms cancel out and we get the expected result.

Numerical scheme for the piano model

The numerical scheme for the piano problem is obtained by "putting together" all the schemes previously described for each part of the instrument, namely (96) (hammer -strings), (105, 109, 113, 114) (stringsoundboard) and (125) (soundboard -air). This gives the following system :

                                                                                                                          m H ξ n+1 h -2ξ n h + ξ n-1 h ∆t 2 = - i F n i,h , F n i,h = -k H δΨ + H d n+1 i,h , d n-1 i,h + r H Φ + H (d n+1 i,h ) -Φ + H (d n-1 i,h ) 2∆t , d k i,h = ξ + ξ k i,h -ξ k h , ξ k i,h = L H,h • Q k i,h , k = n -1, n + 1 M s h Q n+1 i,h -2Q n i,h + Q n-1 i,h ∆t 2 + R s h Q n+1 i,h -Q n-1 i,h 2∆t + A s,D i,h {Q i,h } n θ + A s,p i,h {Q i,h } n 1 12
+ ∇U (Q

n+ 1 2 i,h , Q n-1 2 i,h ) = F n i,h L H,h + F p,n i,h L p,h + F p,n i,h L r,h (130a)
L p,h • Q n+1 i,h -Q n-1 i,h 2∆t = L s,h • U n+ 1 2 p,h -U n-1 2 p,h ∆t , L r,h • Q n+1 i,h -Q n-1 i,h 2∆t = 0 (130b) U n+ 1 2 p,h = U n p,h (t n+ 1 2 ), Un+ 1 2 p,h = d U n p,h dt (t n+ 1 2 ), where (130c)
         d 2 dt 2 U n p,h (t) + Λ p d,h d dt U n p,h (t) + Λ p h U n p,h (t) = - i F p,n i,h . L s,h + C t pp,h P n+ 1 2 h + P n-1 2 h 2 , in t n-1 2 , t n+ 1 2 U n p,h (t n-1 2 ) = U n-1 2 p,h , d dt U n p,h (t n-1 2 ) = Un-1 2 p,h (130d)
         M p h P n+ 1 2 h -P n-1 2 h ∆t + C h V n h = -C pp,h U n+ 1 2 p,h -U n-1 2 p,h ∆t M V h V n+1 h -V n h ∆t -C t h P n+ 1 2 h = 0 (130e)
The proof of existence and uniqueness of the solution to this discrete scheme has not been established. It faces the same difficulties as the ones emphasized in remark 7. However, we are able to state a general stability result.

Proposition 6. Any solution to (130) satisfies the following energy identity:

E n+ 1 2 h,s,p,a,h -E n-1 2 h,s,p,a,h ∆t = - i k H Φ + H (d n+1 i,h ) -Φ + H (d n-1 i,h ) 2∆t d n+1 i,h -d n-1 i,h 2∆t
- i R s h Q n+1 i,h -Q n-1 i,h 2∆t • Q n+1 i,h -Q n-1 i,h 2∆t - t n+ 1 2 t n-1 2 Λ p d,h d dt Ũn p,h (s) • d dt Ũn p,h (s) ds (131)
where the global energy for the piano is defined by

E n+ 1 2 h,s,p,a,h = E n+ 1 2 h,s,h + E n+ 1 2 p,h + E n+ 1 2 a,h (132)
where the energy E by (121). As a consequence, if θ ≥ 1/4 and under the stability conditions (79) and (122), the scheme (130) is stable.

We shall not detail the proof below, which is quite similar to proofs given in previous sections.

Numerical, algorithmic and computational aspects

In this section we provide some details about the algorithmic aspects and the implementation of the simulation code for the piano based on the numerical schemes of section (4.6). The size of the problem implies some choices in order to ensure the feasibility and performances in terms of computational time and memory storage. In particular, we have designed an algorithm that allows at each time step to update independently the unknowns associated with each physical subsystems (hammer-strings, soundboard, acoustics), the coupling being ensured by solving auxiliary problem of small size with respect to the total size of the problems. Thanks to the linearity of the bridge -soundboard -air propagation chain, this is achieved via a Schur complement procedure which results in the algorithm described in 5.1 whose justification is given in 5.2.

Description of the algorithm

Initialization phase.

Prior to any time iteration, we shall compute the constant coefficients and matrices

           F s = -R ∆t L s,h (133a)
S p = I + 1 2 R ∆t C t pp,h (M p h ) -1 C pp,h (133b)
β p,h = L s,h • S -1 p F s (133c) Time iterations.
Assuming that all discrete unknowns are known up to time t n , the vectors

ξ n+1 h , Q n+1 i,h , U n+ 1 2 p,h , Un+ 1 2 p,h , P n+ 1 2 h , V n+1 h are computed in 5 steps.
Step 1. Compute the auxiliary coefficient α

n-1 2 p,h : α n-1 2 p,h = L s,h • S -1 p F n-1 2 a + F n-1 2 p -U n-1 2 p,h (134)
where

     F n-1 2 a = 1 2 R ∆t C t pp,h 2 P n-1 2 h -∆t (M p h ) -1 C h V n h F n-1 2 p = S 0 ∆t U n-1 2 p,h + S 1 ∆t Un-1 2 p,h + 1 2 R ∆t C t pp,h (M p h ) -1 C pp,h U n-1 2 p,h (135)
Step 2. Update the { hammer -strings } unknowns ξ n+1 h , Q n+1 i,h , together with the forces at the bridge F p,n i,h and F r,n i,h , by solving the nonlinear system:

                                                 m H ξ n+1 h -2ξ n h + ξ n-1 h ∆t 2 = - i F n i,h , F n i,h = -k H δΨ + H d n+1 i,h , d n-1 i,h + r H Φ + H (d n+1 i,h ) -Φ + H (d n-1 i,h ) 2∆t , (136a)
d k i,h = ξ + ξ k i,h -ξ k h , ξ k i,h = L H,h • Q k i,h , k = n -1, n + 1 (136b) M s h Q n+1 i,h -2Q n i,h + Q n-1 i,h ∆t 2 + R s h Q n+1 i,h -Q n-1 i,h 2∆t + A s,D i,h {Q i,h } n θ + A s,p i,h {Q i,h } n 1 12 (136c) + ∇U (Q n+ 1 2 i,h , Q n-1 2 i,h ) = F n i,h L H,h + F p,n i,h L p,h + F p,n i,h L r,h (136d)
L p,h • Q n+1 i,h -Q n-1 i,h 2∆t = 1 ∆t α n-1 2 p,h + i F p,n i,h β p,h , L r,h • Q n+1 i,h -Q n-1 i,h 2∆t = 0 (136e)
Here the difficulty lies in the nonlinear nature of this system. The Jacobian matrix of this nonlinear system is expressed analytically [START_REF] Chabassier | Modélisation et simulation numérique d'un piano par modèles physiques[END_REF] and is used in a "modified" Newton -Raphson method, for which the convergence properties are strongly linked to the choice of a guess point, and the choice of a stopping criterion (see [START_REF] Bertsekas | Nonlinear Programming: 2nd Edition[END_REF]). The modification consists of updating the Jacobian matrix of the standard Newton method only when the number of iterations becomes too large.

Step 3. Update the value of U n+ 1 2 p,h by solving

S p U n+ 1 2 p,h = F n-1 2 a + F n-1 2 p + i F p,n i,h F s (137)
where

F n-1 2 a , F n-1 2 p
have been defined in (135), F s in (133a). Notice that this matrix is dense and that its size is equal to the number of considered modes. Since several inversions are required at each time step, the Cholesky factorization of the symmetric matrix R

-1 2 ∆t S p R 1 2
∆t is done once and for all before all the time iterations.

Step 4. Update the acoustic unknowns:

     P n+ 1 2 h = P n-1 2 h -∆t (M p h ) -1 C h V n h -(M p h ) -1 C pp,h U n+ 1 2 p,h -U n-1 2 p,h V n+1 h = V n h + ∆t M V h -1 C t h P n-1 2 h (138)
This is an explicit step thanks to the diagonal nature of mass matrices, but is however the most costly one due to the size of the computational 3D domain. This is overcome via a parallel implementation related to a geometrical partitioning of the unknowns.

Step 5. Update the value of Un+ 1 2 p,h using the analytic formula:

Un+ 1 2 p,h = Ṡ0 ∆t U n-1 2 p,h + Ṡ1 ∆t Un-1 2 p,h + Ṙ∆t - i F p,n i,h L s,h + C t pp,h P n+ 1 2 h + P n-1 2 h 2 (139)
Notice that, if a few modes are considered on the plate (M is small), the computational cost is dominated by step 4 which mainly corresponds to the resolution of standard acoustic equations with a source term concentrated on the plate.

Justification of the algorithm

Step 4 and 5 of the previous algorithm are obvious. To obtain step 3 we use the analytic formula (106) in order to express U n+ 1 2 p,h as a function of the soundboard's past and the right hand side:

U n+ 1 2 p,h = S 0 ∆t U n-1 2 p,h + S 1 ∆t Un-1 2 p,h + R ∆t - i F p,n i,h L s,h + C t pp,h P n+ 1 2 h + P n-1 2 h 2 (140)
This formula does not define U

n+ 1 2 p,h explicitly since P n+ 1 2 h
is unknown. This is why we use (130e) to obtain:

P n+ 1 2 h = P n-1 2 h -∆t (M p h ) -1 C h V n h -(M p h ) -1 C pp,h U n+ 1 2 p,h -U n-1 2 p,h . (141)
This formula is used to obtain:

I + 1 2 R ∆t C t pp,h (M p h ) -1 C pp,h Sp U n+ 1 2 p,h = F n-1 2 a + F n-1 2 p + i F p,n i,h F s (142)
where

           F n-1 2 a = 1 2 R ∆t C t pp,h 2 P n-1 2 h -∆t (M p h ) -1 C h V n h F n-1 2 p = S 0 ∆t U n-1 2 p,h + S 1 ∆t Un-1 2 p,h + 1 2 R ∆t C t pp,h (M p h ) -1 C pp,h U n-1 2 p,h F s = -R ∆t L s,h
which gives equation (137) of step 3.

Formula (142) defines U n+ 1 2 p,h as an affine function of F p,n i,h , which is unknown. Let us exploit this feature in the strings -soundboard first coupling conditions (130b):

L p,h • Q n+1 i,h -Q n-1 i,h 2∆t = 1 ∆t L s,h • S -1 p F n-1 2 a + F n-1 2 p -U n-1 2 p,h α n-1 2 p,h + i F p,n i,h 1 ∆t L s,h • S -1 p F s β p,h (144)
where the scalar number β p,h does not change as time goes by while the scalar number α n-1 2 p,h has to be updated at each time step. This is a closing condition for the { hammer -strings -forces at the bridge } system. This justifies steps 1 and 2.

Numerical results

In this section we present various numerical results of a simulation of a realistic piano which corresponds to a grand Steinway D, for which physical measurements of the first part of this paper [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF] have been performed. We refer the reader to the technical report [START_REF] Chabassier | Physical parameters for piano modeling[END_REF] for the numerical values corresponding to this piano, of all the physical coefficients involved in the model.

These results are split into three subsections. First, we illustrate how the code can be used to visualize the main physical phenomena involved in a piano sound generation. Next, we compare simulations obtained with different string models, that aim at justifying our final choice. Finally, the numerical results are compared to experimental data. Simulated sound pressures will be evaluated at a point close to the piano, vertical to point P of figure 4. P (0.85, 1.459) (0,0) In our presentation we shall emphasize the fact that our model is able to fulfill the five objectives we aimed at in the first part of this paper [CCJ14], except for objective 3 which has been achieved by construction of our model. Let us recall them for ease of reading.

Objective 1. Represent attenuation phenomena which are selective in frequency.

Objective 2. Integrate non linearities in order to discriminate piano and fortissimo sounds.

Objective 3. Account for some mechanism of transmission of the longitudinal string's displacements to the soundboard.

Objective 4. Reproduce the inharmonicity effects.

Objective 5. Account for the phantom partials.

Numerical illustration

The following simulations have been obtained for note C2, with an initial velocity of 4 m/s corresponding to a forte dynamics. This note corresponds to a triplet of 1,6 meter long strings (N s = 3).

The physical computational domain is obtained by extruding the 2d quadrilateral mesh displayed in figure 3(a) by 16 layers, which gives an hexahedral mesh for the parallelepipedic box of volume 2.2 × 3.2 × 0.4 m 3 . The rim is 30 centimeter high and is a whole in the mesh. To simulate an infinite exterior domain, we use perfectly matched layers (which increases slightly the computational domain) [START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF]. We use fourth order polynomials for all Galerkin approximations (r a = r p = r s = 4). One degree of freedom is put every 2 millimeters on the string. This gives about 800 geometrical degrees of freedom on the string (namely 2,400 unknowns). We use 141,329 geometrical degrees of freedom on the soundboard (i.e. 423,987 unknowns) and 269,532 hexahedra in the acoustical domain (without PML layers), which gives about 17 millions of unknowns for the pressure and about 101 millions for the acoustical velocity. The number of computed modes on the soundboard is M =2,400 leading to a maximal eigenfrequency of about 10 kHz. Concerning the time discretisation parameters, we take θ = 1/4, and the time step is ∆t = 10 -6 second. With this choice, stability conditions (122) and (79) are satisfied. Simulations have been run on parallel platforms. 24 hours of real time on 300 processors are needed to provide the piano sound generation during one second. In figure 5 is displayed the time evolution, during the 50 first milliseconds after hammer blow, of the energies of each subsystem and of the hammer strings interaction force. The strings interact with the hammer only in the 5 first milliseconds of the simulation, and we can observe the progressive transfer of the energy from the strings to the acoustic region through the soundboard. In figure 6 are displayed on the 3d pictures the evolution of the acoustic pressure in two vertical slices, and of the soundboard elevation in the horizontal plane. Above each 3d picture, the string deformation is displayed and the longitudinal displacement is represented in color scale. Notice that the longitudinal wave propagates faster than the transversal one, hence the deformation of the soundboard begins before the string deformation reaches the attach point. This illustrates the presence of a precursor, which can in fact be attributed to the nonlinear behavior of the string, see section 6.2. Finally, we have computed two sounds corresponding to note C♯5 for two different initial hammer velocities : 0.5 m/s (piano dynamics) and 4 m/s (fortissimo dynamics). In figure 7 are displayed the spectrograms of the associated simulated sound pressure 18 cm above point P of figure 4, which are clearly different and emphasize the nonlinear behavior of the instrument. Moreover, high frequencies are more rapidly attenuated than low frequencies. This shows that we reached our objectives 1 and 2.

Influence of the string model

Our goal is to analyze the influence of strings' model on the resulted computed sound. More precisely, we shall compare:

• Model 1: the stiff linear model in (u i , ϕ i ) of section (3.2) of [CCJ14],

• Model 2: the geometrically exact model in (u i , v i ) of section (3.3) of [CCJ14],

• Model 3: the stiff nonlinear model in (u i , v i , ϕ i ), presented in this paper.

In figure 8 we first compare these models for note C2, through the spectra of the computed pressure field 1 cm above point P of figure 4, in the frequency range [6, 9] kHz. Thick peaks can be observed in all three figures, which can be associated with the eigenfrequencies of the corresponding linearized models. As expected, these are the same for models 3 and 1 (represented with solid lines): one recovers the inharmonic spectrum for the unknown u i corresponding to the eigenfrequencies (f n ) n≥1 described in section 3.2.1 of [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF], see also equation (71). For model 2, one recovers (in dashed lines) the spectrum (f Ale n ) n≥1 of d'Alembert's model described in section 3.3.3 of [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF], which differs from the inharmonic spectra of models 1 and 3 more and more as the frequency is increased. In addition, additional peaks (indicated by magenta circles) correspond to the linear vibrations in the longitudinal direction for the unknown v i and coincide with the eigenfrequencies (f ℓ n ) n≥1 described in section 3.2.1 of [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF], see also equation (71). These eigenfrequencies are much fewer

Figure 1 .Figure 2 .

 12 Figure 1. Geometric configuration of the piano.

  (a) 2D mesh to be extruded delimiting the soundboard, the rim, and the acoustic region (b) Trace of the volumic mesh on the obstacles to sound propagation : the soundboard and the rim.

Figure 3 .

 3 Figure 3. Meshes.

n+ 1 2 h

 2 ,s,h is defined by (99), the energy E n+ 1 2 p,h by (110) and the energy E

Figure 4 .

 4 Figure 4. Simulated sound pressures are evaluated at a point vertical to point P.

  Energy vs time for note C2 (three strings). Black solid line: Total energy. Blue --line: Hammer. Green -•line: Strings. Red • • • line: Soundboard. Turquoise thick -line: Air. Left: Linear scale. Right: Logarithmic scale.

Figure 5 .

 5 Figure 5. Note C2, energy transfers and hammer string interaction force.

Figure 6 .

 6 Figure 6. Time evolution of some variables of the piano model for string C2. The transverse displacement of the string is represented in the upper parts of the figures, while the longitudinal displacement is shown through shading in the string thickness (upper scale). The displacement of the soundboard is shown in the lower parts, while the pressure is shown in two vertical planes which cross at the point where the string is attached to the bridge: x 1 = 0.59 m and x 2 = 1.26 m. The lower scale is related to the sound pressure. The scale of the soundboard's displacement is adjusted over time in order to see the evolution of waves clearly. (a) t = 0.4 ms. (b) t = 1.1 ms. (c) t = 2.1 ms. (d) t = 3.1 ms. (e) t = 4.1 ms. (f) t = 5.1 ms. (g) t = 7.1 ms. (h) t = 8.1 ms. (i) t = 16.1 ms.

  C♯5 note, p dynamics (hammer velocity 0.5 m/s) C♯5 note, ff dynamics (hammer velocity 4.5 m/s)

Figure 7 .

 7 Figure 7. Spectrogram (time -frequency diagram) of simulated note C♯5 in piano (Left) and fortissimo (Right) dynamics.

Figure 8 .Figure 9 .

 89 Figure 8. Power spectrum of simulated sound pressures using different string models.

  Spectrogram of the simulated sound pressure

Figure 10 .

 10 Figure 10. Note C♯5 : comparison between recorded and simulated sounds. The cyan arrow points out the longitudinal partial.

Figure 11 .

 11 Figure 11. Measured (left) and simulated (right) starting transients of hammer acceleration, string transversal displacement and sound pressure for note D♯1.

because longitudinal waves are much faster.

With model 3, additional secondary peaks (of lower amplitude) are present, indicated by orange triangles. These are the fantom partials referred to in the literature [START_REF] Conklin | Piano strings and "phantom" partials[END_REF], see also section 2 of [START_REF] Chabassier | Time domain simulation of a piano. Part I : Model description[END_REF]. Their frequencies correspond to sums or differences of two different frequencies of the inharmonic spectrum (f n ) n≥1 , and can be attributed to the quadratic nonlinearities of model 3. In spite of the nonlinearity of model 2, such fantom partials cannot be observed because its linearized spectrum (f Ale n ) n≥1 is harmonic. This demonstrates that only model 3 is able to reach both objectives 4 and 5.

In figure 9 are displayed the first milliseconds of sound (pressure field) for each string model 1, 2 and 3. The sound precursors can be clearly observed on the last two curves that are linked to the nonlinear nature of the model and the longitudinal vibrations (the unknown v i ). This illustrates the fact that the bridge condition is able to model the transmission of longitudinal waves to the soundboard, which was the object of our objective 3.

Comparison between measurements and simulations

In this last paragraph, we present comparisons between simulated and recorded signals. The experiments have been carried out on a grand Steinway D provided by IRCAM. In figure 10, the computed and recorded waveforms for note C♯5 : the pressure fields evaluated 18 cm above point P of figure 4 with respect to time. Despite the fact that the chosen points are not exactly the same, which is a source of discrepancies between the two curves, one observes a very similar behavior, in particular the double decay phenomenon(see [START_REF] Weinreich | Coupled piano strings[END_REF]). There is also a very good agreement between the two associated seismograms, in particular, in the attenuation of the partials. Notice the longitudinal partial present in both spectrograms and indicated with a cyan arrow. They are distinguishable from the main series of partials for two principal reasons : first, their attenuation follows a different law with respect to the frequency, which makes them longer than their neighbors. This is due to the fact that the longitudinal wave is transmitted less efficiently to the soundboard than the transversal wave. Second, they are thinner since the longitudinal frequencies only depend on the length, the Young's modulus, and the density of the string, and not on the tension at rest which is the only difference between the three strings of the choir in our model (and leads to three slightly different series of transversal partials). The value of the longitudinal partial is not exactly the same in the recorded and simulated sounds, because the strings parameters of the simulation have not been precisely adjusted.

Finally, in figure 11 are displayed the waveforms for note D♯1 of the hammer acceleration, the transversal string displacement at a point located 1.749 m from the agraffe, and the pressure field 18 cm above point P of figure 4. We observe a very good qualitative agreement between simulated (right) and recorded (left) curves. The differences between hammer accelerations are due to the fact that the hammer shank has not been taken into account in our model (see [START_REF] Chabassier | Energy based simulation of a Timoshenko beam in non-forced rotation. Influence of the piano hammer shank flexibility on the sound[END_REF] for more details). It is to our knowledge the first piano simulations to reach this quality of restitution.

Acknowledgments

Experiments presented in this paper were carried out using the PLAFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr/). Computer time for this study was also provided by the computing facilities MCIA (Mésocentre de Calcul Intensif Aquitain) of the Université de Bordeaux and of the Université de Pau et des Pays de l'Adour.