Juliette Chabassier

Antoine Chaigne

Patrick Joly

Modeling and simulation of a grand piano

Keywords: piano, modeling, time-domain, simulation, nonlinear strings, string-soundboard coupling, radiation piano, modélisation, domaine temporel, simulation, cordes non linéaires, couplage cordes -table d'harmonie, rayonnement

A time-domain global modeling of a grand piano is presented. The string model includes internal losses, stiffness and geometrical nonlinearity. The hammer-string interaction is governed by a nonlinear dissipative compression force. The soundboard is modeled as a dissipative bidimensional orthotropic Reissner-Mindlin plate where the presence of ribs and bridges is treated as local heterogeneities. The coupling between strings and soundboard at the bridge allows the transmission of both transverse and longitudinal waves to the soundboard. The soundboard is coupled to the acoustic field, whereas all other parts of the structure are supposed to be perfectly rigid. The acoustic field is bounded artificially using perfectly matched layers (PML). The discrete form of the equations is based on original energy preserving schemes. Artificial decoupling is achieved, through the use of Schur complements and Lagrange multipliers, so that each variable of the problem can be updated separately at each time step. The capability of the model is highlighted by series of simulations in the low, medium and high register, and through comparisons with waveforms recorded on a Steinway D piano. Its ability to account for phantom partials and precursors, consecutive to string nonlinearity and inharmonicity, is particularly emphasized.

Introduction

In this paper, an extensive and global model of a piano is presented. Its aim is to reproduce the main vibratory and acoustic phenomena involved in the generation of a piano sound from the initial blow of the hammer against the strings to the radiation from soundboard to the air. Compared to previous studies, the prime originality of the work is due to the string model which takes both geometrical nonlinear effects and stiffness into account. Other significant improvements are due to the combined modeling of the three main couplings between the constitutive parts of the instrument: hammer-string, string-soundboard and soundboard-air coupling.

Although a vast literature exists on the piano and its subsystems (strings, hammer, soundboard, radiated field), there are only a few examples of a complete computational model of the instrument. One noticeable exception is the work by Giordano and Jiang [START_REF] Giordano | Physical modeling of the piano[END_REF] describing the modeling of a linear string coupled to soundboard and air using finite differences. Compared to this reference, our work is based on a more accurate description of the piano physics, and also pays more attention to the properties of the numerical schemes used for solving the complex system of coupled equations. The strategy used here for the piano is similar to the one developed previously by two of the authors for the guitar [START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF] and the timpani [START_REF] Rhaouti | Time-domain modeling and numerical simulation of a kettledrum[END_REF]. The physical model is composed of a set of equations governing the hammer-string contact and the wave propagation in strings, soundboard and air, in the time domain. The input parameters of the equations are linked to the geometry and material properties of the propagation media. The equations are then discretized in time and space, in order to allow the numerical resolution of the complete system. Original numerical schemes are developed in order to ensure stability, sufficient accuracy and the conservation of energy. In this respect, this strategy has a direct continuity with the work by Bilbao [START_REF] Bilbao | Conservative numerical methods for nonlinear strings[END_REF].

The numerical dispersion of the schemes is maintained sufficiently small so that the difference between real and simulated frequencies are comparable to the just noticeable differences of the human ear in the audible range (around 1 %). The validity of the numerical model is assessed by careful analysis, in time and frequency, of the most significant variables of the problem: hammer force, string, bridge and soundboard displacements (resp. velocities or accelerations), sound pressure. In a first step, using typical values of parameters found in the literature, the numerical results are expected to reproduce the main features of piano sounds, at least qualitatively. In a second step, the input parameters are the results of measurements on real pianos, and a more thorough comparison is made between real and simulated tones.

One motivation at the origin of this study was to reproduce the effects of the geometrical nonlinearity of piano strings both in time (precursors [START_REF] Askenfelt | Observations on the transient components of the piano tone[END_REF], transverse-longitudinal coupling) and frequency domain (phantom partials [START_REF] Conklin | Generation of partials due to nonlinear mixing in a stringed instrument[END_REF]). The purpose of the simulations is to get a better understanding of these phenomena: the aim is to exhibit quantitative links between some control parameters (initial hammer velocity, string tension,...) and observed waveforms and spectra. Experimental observations of string and soundboard spectra recorded on real pianos suggest, in addition, that the string-soundboard coupling plays a crucial role in the transmission of longitudinal string components [START_REF] Giordano | Motion of a piano string: Longitudinal vibrations and the role of the bridge[END_REF], and thus attempting to reproduce and understand these effects is another challenge.

In pianos, it is well-known that the initial transients are perceptually highly significant. One major attribute of spectral content during the attack is due to the presence of soundboard modes excited by the string force at the bridge. The low-frequency modes, in particular, have low damping: they are well separated and clearly recognizable on the spectra [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part II. Piano structure[END_REF]. An accurate model of stringsoundboard coupling has the capability of accounting for such effects. Soundboard modeling also allows to explore the effects of bridge and ribs distribution on the produced sound.

Finally, air-soundboard coupling is a necessary step for simulating sound pressure. Because of the wideband modeling, the computation of a 3-D pressure field is highly demanding. High performance parallel computing was necessary here. However, the results are very valuable, since one can make direct auditory comparisons with real piano sounds. Simulation of a 3-D field also yields information on directivity and sound power.

The physical model of the piano used in this work is presented in Section 2. Emphasis is put on what we believe to be the most innovative parts of the model: the nonlinear stiff string and its coupling with an heterogeneous orthotropic soundboard via the bridge. Hammer-string interaction and air-soundboard coupling are also described. In Section 3, the general method used for putting the model into a discrete form is explained. The properties of the main numerical schemes retained for each constitutive part of the instrument and the coupling conditions are given without demonstrations, with explicit references to other papers published by the authors, that are more oriented on the mathematical and numerical aspects of the model.

The validity of the model is evaluated in Section 4 through analysis of a selection of simulated piano notes in the low, medium and treble register, respectively. The motion of the main elements (piano hammer, strings, bridge, soundboard, air) are analyzed alternatively in time, space, and/or frequency domain. Because of the important effort put on the nonlinear string, the capability of the model to reproduce amplitude-dependent phenomena is emphasized. Comparisons are also made with measurements performed on a Steinway grand piano (D model). In these series of experiments, the motion of hammer, string, bridge and soundboard, and the sound pressure were recorded simultaneously. This shed useful light on the transmission and transformation of the signals from hammer to sound and, through comparisons, on the ability of the model to account for the coupling between the constitutive parts of the piano.

Presentation of the model

During the development of the model, particular attention was paid to the following guidelines:

• The prime function of the piano is to generate sounds. However, its conception is not only governed by dynamical and acoustical laws, but also by purely static and mechanical considerations, as it is the case, for example, for the design of frame, rim and feet. In addition, from the point of view of the player, the piano must fulfill a number of playability rules. Our purpose is not to mimic piece by piece the construction of a real instrument, but rather to describe the major phenomena that are at the origin of piano tones. This is achieved here by considering a set of equations governing the nonlinear hammer-string contact, the wave propagation in strings, the vibrations of the heterogeneous orthotropic soundboard radiating in air, and the reciprocal coupling between strings and soundboard at the bridge.

• Several simplifying assumptions were made in the model. The hammer is supposed to be perfectly aligned with the strings. The agraffe is assumed to be rigidly fixed with a simply supported end condition for the strings. Both the string-soundboard and soundboard-air couplings are supposed to be lossless. The soundboard is considered as simply supported along its edge, and the "listening" room is anechoic with no obstacle except the piano itself. Finally, the action of the mechanism prior to the shock of the hammer against the strings is ignored. In other words, the tone starts when the hammer hits the strings with an imposed velocity.

• Most of the physical parameters of hammers, strings and soundboards included in the model are obtained from standard string scaling and geometrical data from piano, strings and soundboards manufacturers. When necessary, these data were complemented by results of our own measurements on various instruments. However, some difficulties arose for the physical description of damping models in hammers, strings and soundboards since the underlying physical mechanisms are often very complex, and because an accurate analytical description of these phenomena are not available in the present state-of-the-art. Thus, it was decided to use here approximate models based on experimental data, as it is commonly done in other fields of structural dynamics [START_REF] Balmes | Modeling damping at the material and structure level[END_REF].

• As shown in Section 3, the numerical formulation of the model is based on a discrete formulation of the global energy of the system, which ensures stability. This requires that the continuous energy of the problem is decaying with time. The global model of the piano is thus written according to this condition.

In what follows, an energy decaying model is developed for each constitutive part of the instrument, successively. Finally, a global piano model is proposed.

Strings

The selected string model accounts for large deformations, inducing geometrical nonlinearities, and intrinsic stiffness. These are two essential physical phenomena in piano strings [START_REF] Jarvelainen | Audibility of the timbral effects of inharmonicity in stringed instrument tones[END_REF][START_REF] Bank | Perception of longitudinal components in piano string vibrations[END_REF]. The governing equations correspond to those of a Timoshenko beam under axial tension. Modeling the stiffness with a Timoshenko model, rather than with an Euler-Bernoulli model as in a previous study by Chaigne and Askenfelt [START_REF] Chaigne | Numerical simulation of piano strings. I. A physical model for a struck string using finite-difference methods[END_REF], is motivated by both physical and numerical reasons. As shown in Fig. 13 of Section 4, the physical dispersion predicted by the Timoshenko model shows a good agreement with experimental data in a wider frequency range than the Euler-Bernoulli model. It also yields an asymptotic value for the transverse wave velocity as the frequency increases, in contrast with the Euler-Bernoulli model. This latter property is not only more satisfying from the point of view of the physics, but it is also more tractable in the simulations. The geometrical nonlinearities are described using a standard model [START_REF] Morse | Theoretical Acoustics[END_REF].

Let us call ρ the density of the string, A the area of its cross-section, T 0 its tension at rest, E its Young's modulus, I its stiffness inertia coefficient, G its shear coefficient, and κ the nondimensional Timoshenko parameter. We denote u s (x, t) the transverse vertical displacement, v s (x, t) the longitudinal displacement, and ϕ s (x, t) the angle of the cross-sections with the plane normal to the string (see Fig. 1). The space variable is x ∈ [0, L], where L is the length of the string at rest. For the end conditions, we assume zero displacement (in both transverse and longitudinal directions) and zero moment. These conditions are motivated by usual observations at the agraffe (x = 0) [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part III. Piano strings and scale design[END_REF], and will be revisited later at point x = L when considering the coupling with the soundboard. Finally, the string is considered at rest at the origin of time. We get the following system [START_REF] Giordano | Physical modeling of the piano[END_REF], where S is a source term. As shown in the next paragraph, this term is supposed to account for the action of the hammer against the strings. The model is written:

String motion:

(1)

                                           ρA ∂ 2 u s ∂t 2 - ∂ ∂x       EA ∂u s ∂x - (EA -T 0 ) ∂u s ∂x ∂u s ∂x 2 + 1 + ∂v s ∂x 2       + AGκ ∂ ∂x ϕ s - ∂u s ∂x = S, ρA ∂ 2 v s ∂t 2 - ∂ ∂x       EA ∂v s ∂x - (EA -T 0 ) 1 + ∂v s ∂x ∂u s ∂x 2 + 1 + ∂v s ∂x 2       +AGκ ∂ ∂x ϕ s - ∂u s ∂x = 0, ρI ∂ 2 ϕ s ∂t 2 -EI ∂ 2 ϕ s ∂x 2 + AGκ ∂ ∂x ϕ s - ∂u s ∂x = 0
Boundary conditions:

(2) A distinction should be made in the system of Eqs. (1) between plain and wrapped strings. For plain strings, the density ρ, cross-sectional area A, and stiffness coefficient I refer to the parameters of the metallic wire (usually steel), without ambiguity. For wrapped strings, A = A c is the cross-sectional area of the core, I = I c is the stiffness coefficient of the core, whereas ρ = ρ c F where ρ c is the density of the core, and F is a wrapping factor (see Conklin [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part III. Piano strings and scale design[END_REF]) defined by the expression:

u s (x = x b , t) = v s (x = x b , t) = ∂ϕ s ∂x (x = x b , t) = 0, for x b ∈ {0, L}, Initial conditions: (3)      u s (x, t = 0) = v s (x, t = 0) = ϕ s (x, t = 0) = 0, ∂u s ∂t (x, t = 0) = ∂v s ∂t (x, t = 0) = ∂ϕ s ∂t (x, t = 0) = 0.
F = 1 + ρ w A w ρ c A c , (4)
where A w is the cross-sectional area of the wrapping with density ρ w (usually copper). In both cases, the other parameters E, T 0 , G, κ refer to the material of the wire (resp. core). In short, the wrapping affects the inertial forces only.

The eigenfrequencies of the linearized version of system (1) associated with boundary conditions (2) can be computed analytically (see, for example [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]), which yield three different series. The following expressions hold for partial's rank ℓ ≥ 1. The flexural eigenfrequencies are given by:

f trans ℓ = ℓf trans 0 1 + ǫ ℓ 2 + O(ℓ 5 ),
where

f trans 0 = 1 2L T 0 ρA , ǫ = π 2 2L 2 EI T 0 1 - T 0 EA . (5)
The shear eigenfrequencies satisfy:

f shear ℓ = f shear 0 1 + η ℓ 2 + O(ℓ 4 ),
where

f shear 0 = 1 2π AGκ ρI , η = π 2 2L 2 EI + IGκ AGκ . (6)
These shear frequencies, which are situated well above the audio range, are not perceived by the human ear. Finally, the longitudinal eigenfrequencies read

f longi ℓ = ℓf longi 0 , where f longi 0 = 1 2L E ρ . (7)
Eq. [START_REF] Askenfelt | Observations on the transient components of the piano tone[END_REF] shows that the model accounts for inharmonicity. The inharmonicity coefficient ǫ is close to, but slightly different from, the usual coefficient obtained with the Euler-Bernoulli model [START_REF] Fletcher | Normal vibration frequencies of a stiff piano string[END_REF]. The speed of waves associated to equations ( 5) and ( 7) are the transverse speed c trans = T 0 /ρA and the longitudinal speed c longi = E/ρ. For the string C2, for example, we obtain c trans = 209 m•s -1 and c longi = 2914 m•s -1 . Since both types of waves are present on the string (because of nonlinear transverse to longitudinal coupling) near the hammer contact point, the speed ratio explains why the longitudinal waves arrive first at the bridge, inducing a precursor.

With the objective to simulate realistic tones, it is essential to account for the observed frequencydependent damping on strings. Damping phenomena in structures in general, and in strings in particular, are hard to apprehend, for many reasons (lack of available measurements, uneasy separation of the mechanisms, difficulty in modeling microscale phenomena such as dislocations, in the equation of motion,. . . ). Therefore, following the strategy used in previous studies, we use a simple model that globally and approximately accounts for the damping effect, without pretending to model the underlying physics in details [START_REF] Bensa | The simulation of piano string vibration: From physical models to finite difference schemes and digital waveguides[END_REF]. In practice, viscoelastic terms are added under the form:

2 ρA R u ∂u s ∂t -2 T 0 η u ∂ 3 u s ∂t ∂x 2 (8)
where R u and η u are empirical damping coefficients to be determined from measured sounds for each string, through comparisons between simulated and measured spectrograms. The resulting damping law in the frequency domain is the sum of a constant term 2 ρA R u and a quadratic term 2 T 0 η u f 2 ℓ . Similar damping terms have to be added for shear and longitudinal motions in order to avoid unpleasant endless frequencies in the simulated tones.

2 ρA R v ∂v s ∂t -2 EA η v ∂ 3 v s ∂t ∂x 2 , 2 ρI R ϕ ∂ϕ s ∂t -2 EI η ϕ ∂ 3 ϕ s ∂t ∂x 2 . (9)
Again, the order of magnitude for (R v , R ϕ ) and (η v , η ϕ ) are determined from experiments and/or trial-and-error procedures, considering the absence of any other reliable method.

When the system is subjected to a source term S(x, t) in the transverse direction, it is possible to show the following energy identity [START_REF] Chabassier | Time domain simulation of a piano. Part I: model description[END_REF]:

dE s dt ≤ L 0 S ∂u s ∂t dx where E s (t) = E s,kin (t) + E s,pot (t) (10)
E s,kin (t) = ρA 2 L 0 ∂u s ∂t 2 dx + ρA 2 L 0 ∂v s ∂t 2 dx + ρI 2 L 0 ∂ϕ s ∂t 2 dx E s,pot (t) = T 0 2 L 0 ∂u s ∂x 2 dx + EA 2 L 0 ∂v s ∂x 2 dx + EI 2 L 0 ∂ϕ s ∂x 2 dx + AGκ 2 L 0 ∂u s ∂x -ϕ s 2 dx + (EA -T 0 ) L 0   1 2 ∂u s ∂ x 2 + 1 + ∂v s ∂x - ∂u s ∂x 2 + 1 + ∂v s ∂x 2   dx
This quantity is always positive in practice since, for real piano strings, we have EA > T 0 . Advantage of this property will be taken in Section 3 to derive stable numerical schemes for solving the nonlinear system of equations.

Hammer

We now turn to the interaction between the hammer and the strings. Depending on the note, the hammer strikes one, two or three strings. The components of the motion of the ith string's of a given note's set are written (u s,i , v s,i , ϕ s,i ). Since the strings belonging to the same note are slightly detuned [START_REF] Weinreich | Coupled piano strings[END_REF], each string has a different tension at rest T 0,i in system [START_REF] Giordano | Physical modeling of the piano[END_REF]. The hammer's center of gravity is supposed to be moving along a straight line orthogonal to the strings at rest. Its displacement on this line is represented by ξ(t). The interaction force between the hammer and the ith-string of a note is distributed on a small portion of the string, through a spreading function δ H (represented in Fig. 2) localized around the impact point x 0 , with

L 0 δ H (x) dx = 1.
As a consequence, we denote u s,i the weighted average of the transverse displacement of the string:

u s,i (t) = L 0 u s,i (x, t) δ H (x) dx. ( 11
)
RR n°8181 The interaction force depends on the distance d(t) between the hammer and the string: if d(t) > ξ, there is no contact and the force is zero. If d(t) ≤ ξ, the force is a function of the distance.

According to previous studies, we define the function:

Φ(d) = ξ -d + p (12)
where (•) + means "positive part of", and where p is a real positive nonlinear exponent. In practice, this coefficient varies between 1.5 and 3.5 [START_REF] Stulov | Dynamic behavior and mechanical features of wool felt[END_REF][START_REF] Giordano | Piano hammers and their force compression characteristics: does a power law make sense ?[END_REF][START_REF] Chaigne | Numerical simulation of piano strings. I. A physical model for a struck string using finite-difference methods[END_REF]. In order to account further for the observed hysteretic behavior of the felt [START_REF] Stulov | Dynamic behavior and mechanical features of wool felt[END_REF], a dissipative term is added in the expression of the force. In summary, the parameters characterizing the mechanical behavior of the hammer are its equivalent mass M H , its stiffness coefficient K H i and its dissipation coefficient R H i . The index i here indicates that the model can eventually account for the fact that the interaction with the hammer might not be uniform for each string of a note. The interacting force between the hammer and the ith-string finally reads

F H i (t) = K H i Φ u s,i (t) -ξ(t) + R H i d dt Φ u s,i (t) -ξ(t) . ( 13
)
The hammer is submitted to the sum of these forces. Conversely, the right hand side S(x, t) of the strings system (1) is replaced with F H i (t) δ H (x). We consider an initial position -ξ 0 and an initial velocity v H 0 for the hammer, while the strings are considered at rest at the origin of time. Defining further Ψ(d) = +∞ d Φ(s) ds, one can show (see Chabassier et al. [START_REF] Chabassier | Time domain simulation of a piano. Part I: model description[END_REF]) that any regular solution to the resulting hammer -strings system satisfies the following energy decay:

dE s,h dt ≤ 0, with E s,h (t) = i E s,i (t) + K H i Ψ u s,i (t) -ξ(t) + M H 2 dξ dt (t) 2 ≥ 0 ( 14
)
where E s,i (t) is the energy defined in Eq. ( 10) for the ith string.

Soundboard

The only vibrating element of the piano case considered in the model is the soundboard, all other parts (rim, keybed, lid, iron frame. . . ) being assumed to be perfectly rigid. In view of the small "thickness over other dimensions" ratio of the piano soundboard, a bidimensional Reissner-Mindlin plate model is considered. This model is the bidimensional equivalent to the linear Timoshenko model. It has been preferred here to the Kirchhoff-Love model, because its yields a better estimate

RR n°8181

for the soundboard motion, in the complete audio range. It also has better mathematical properties [START_REF] Chabassier | Time domain simulation of a piano. Part I: model description[END_REF]. The variables of the motion are the vertical transverse displacement u p (x, y, t) at a current point of coordinates (x, y) of the bidimensional plate ω, and two shear angles θ x,p (x, y, t) and θ y,p (x, y, t). These last two variables account for the deviation of the straight segments of the plate from the normal to the medium surface, in the (e x , e y )-referential plane (see Fig. 3). The vector θ p (x, y, t) groups these two angles. The bridge and ribs are considered as heterogeneities of the soundboard, and the orientation of the orthotropy axes can be space dependent. As a consequence, the physical coefficients representing the density ρ p , the thickness δ, the Young's moduli in the two main directions of orthotropy E x and E y , the shear moduli in the three main directions of orthotropy G xy , G xz and G yz , the Poisson's ratios ν xy and ν yx , and the shear coefficient κ x and κ y are functions of space.

Modeling the thickness as a space dependent variable makes it possible to simulate a diaphragmatic soundboard [START_REF] Bilhuber | The influence of the soundboard on piano tone quality[END_REF] (the thickness varies between 6 and 9 mm in the soundboard, between 9 and 35 mm on the ribs, between 29 and 69 mm on the bridge, and between 29 and 95 mm on the crossing areas of ribs and bridge). The soundboard is assumed to be simply supported [START_REF] Mamou-Mani | Numerical simulation of a piano soundboard under downbearing[END_REF] on its edge ∂ω. Finally, a source term is imposed in the transverse vertical direction. The function of this term is to account for both the string's tension at the bridge (see Section 2.4) and the air pressure jump (see Section 2.5). For simplicity, the rotation of the orthotropy axes in the soundboard is ignored in the presentation of the equations, although this flexibility is possible in the model. The Reissner-Mindlin system that governs the motion of transverse displacement and shear angles is written for (x, y) Table 1: Parameters used for the wooden soundboard: Spruce for the table and the ribs, Beech for the bridge.

ρ p E x E y G xy G xz G yz ν xy (kg•m -3 ) (GPa) (GPa) (GPa) (GPa) (GPa)
               ρ p δ 3 12
∂ 2 θ p ∂t 2 -Div δ 3 12 C ε(θ p ) + δ κ 2 • G • (∇u p + θ p ) = 0 (15a) ρ p δ ∂ 2 u p ∂t 2 -div δ κ 2 • G • (∇u p + θ p ) = f (15b) u p = Cε(θ p ) n = 0 on ∂ω (15c) with C ε =      E x 1 -ν xy ν yx - E y ν xy 1 -ν xy ν yx 0 - E x ν yx 1 -ν xy ν yx E y 1 -ν xy ν yx 0 0 0 2G xy        ε xx ε yy ε xy   where the tensor ε is symmetric, G = G xz G yz , κ 2 = κ 2 x κ 2 y (16)
Here, Div is the divergence operator for tensors: Div(τ ) = ∂ j τ i,j , ε is the linearized strain tensor, div is the divergence operator for R 2 vectors and ∇ is the gradient operator. Table 1 gives the parameters used for the soundboard wood (Spruce for the table and the ribs, Beech for the bridge). A prestress term can be added to this model, if necessary, accounting for the static action of the strings on the curved soundboard. Since this action contributes to reduce the initial crown of the soundboard in such a way that it becomes almost flat in normal use, the curvature and the prestress of the soundboard are ignored here.

With regard to the modeling of damping in the soundboard material, a modal approach has been adopted where the modal damping can be adjusted, mode by mode. This method is justified as long as the damping factor is small compared to the eigenfrequency, and is of current use in structural dynamics [START_REF] Hasselman | Modal coupling in lightly damped structures[END_REF]. The modal amplitudes X n (t) of the nth mode associated to the frequency f n are then solution of the second-order uncoupled damped oscillators equations:

d 2 X n dt 2 + α(f n ) dX n dt + (2πf n ) 2 X n = F n ( 17
)
where α is a positive damping function matching experimental data [START_REF] Ege | Vibroacoustics of the piano soundboard:(non)linearity and modal properties in the low-and mid-frequency ranges[END_REF], and F n is the modal contribution of the transversal source term f .

For any positive damping law it is possible to show the following energy identity:

dE p dt ≤ ω f
E p,pot (t) = ω δ 3 12 C ε(θ p ) : ε(θ p ) dx dy + ω δ κ 2 • G |∇u p + θ p | 2 dx dy
Thus, like for the strings-hammer system, we can conclude that the energy of the soundboard decays with time, after extinction of the source.

Strings-soundboard coupling at the bridge

The main purpose of the bridge is to transmit the vibrations from the strings to the soundboard. Conversely, the motion of the soundboard drives the string at the bridge. As highlighted by the spectral content of the precursor signal [START_REF] Podlesak | Dispersion of waves in piano strings[END_REF], it is essential to model the transmission of both transverse and longitudinal waves of the string to the other parts of the structure. The literature is not very broad concerning this part of the instrument, with a few exceptions [START_REF] Giordano | Motion of a piano string: Longitudinal vibrations and the role of the bridge[END_REF][START_REF] Cuenca | Three-dimensional interaction between strings, bridge and soundboard in modern piano's treble range[END_REF]. Also few experimental data are available, and the wave transmission phenomena occurring at the bridge are still not clearly understood. For these reasons, we decided to describe a plausible, though not fully proved, way for transforming the string longitudinal component into a bridge transverse motion.

The method is based on the observation that the strings form a slight angle α with the horizontal plane due to both bridge height and soundboard curvature (see Fig. 4-(a)). We also assume that the bridge moves in the vertical direction only. This, again, seems to be a reasonable assumption since the static tension of the complete set of strings prevents the bridge to have a significant motion in the direction of the strings. At this point, we are aware of the fact that some authors were able to identify and measure such a motion [START_REF] Giordano | Motion of a piano string: Longitudinal vibrations and the role of the bridge[END_REF], but we must admit that we were no able to exhibit similar features in our measurements. It would be probably also necessary in the future to revisit the assumption of ignoring the horizontal bridge motion perpendicular to the strings, that might induce an horizontal component to the string motion. When the hammer strikes the strings, it gives rise to a transversal wave which, in turn, induces a longitudinal wave, because of nonlinear geometrical coupling. The longitudinal wave travels 10 to 20 times faster than the transverse one, and thus it comes first at the bridge (see Section 2). The resulting variation of tension is oriented in the direction of the string at rest. Because of the angle formed by the string with the horizontal plane, this induces a vertical component of the longitudinal force at the bridge, in addition to the transverse force. In our numerical model, the total bridge force is distributed in space in the soundboard by means of a rapidly decreasing regular function χ ω centered on the point where the string is attached on the soundboard (see Fig. 4-(b)). The associated kinematic boundary conditions are the continuity of string and soundboard velocities in the vertical direction ν, and the nullity of the velocity in the horizontal direction τ at this point (see Fig. 4(a)). These conditions are written formally:

   ∂u s,i ∂t (x = L) ∂v s,i ∂t (x = L)    • ν = ω ∂u p ∂t χ ω ,    ∂u s,i ∂t (x = L) ∂v s,i ∂t (x = L)    • τ = 0. ( 19
)
RR n°8181 In addition, the source term in the system of Eqs. ( 15)-(a) to ( 15)-(c) is given by f = -F b (t) χ ω (x, y), where F b (t) is the bridge force associated to the cinematic conditions written in Eq. ( 19). Ignoring the damping terms, and considering only one string for simplicity, this force is written:

F b (t) = cos(α) EA ∂ x u s + AGκ ∂ x u s -ϕ s -(EA -T 0 ) ∂ x u s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 (x = L, t) + sin(α) EA ∂ x v s + (EA -T 0 ) 1 - 1 + ∂ x v s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 ) (x = L, t) (20)
If the magnitude of the transverse motion remains small enough, it becomes justified under some conditions to derive approximate string models using an asymptotic approach [START_REF] Chabassier | Transitoires de piano et non-linéarités des cordes : Mesures et simulations[END_REF]. Such models were used in the past by different authors for analytical [START_REF] Bank | Generation of longitudinal vibrations in piano strings: From physics to sound synthesis[END_REF] or numerical [START_REF] Bilbao | Conservative numerical methods for nonlinear strings[END_REF] purposes. In our case, we will not use this approximate system for the modeling of the piano, and the mathematical reasons for this choice were given explicitely in a previous paper [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations. Application to the vibrating piano string[END_REF]. However, interesting properties can be derived from the approximate expression of the bridge force:

F b (t) ≈ cos(α) T 0 ∂ x u s + AGκ ∂ x u s -ϕ c +(EA -T 0 )∂ x u s ∂ x v s + (EA -T 0 ) (∂ x u s ) 3 2 (x = L, t) + sin(α) EA ∂ x v s + (EA -T 0 ) (∂ x u s ) 2 2 (x = L, t). ( 21
)
From Eq. ( 21), it can be derived that the quadratic and cubic terms in u s generate double and triple combinations of the transverse eigenfrequencies [START_REF] Bank | Generation of longitudinal vibrations in piano strings: From physics to sound synthesis[END_REF]. Combination of longitudinal and transverse eigenfrequencies can also exist potentially, through the product u s v s . All these combinations correspond to the so-called "phantom" partials.

Sound propagation and structural acoustics

We are now interested in the propagation of piano sounds in free space. The rim (occupying the space Ω r ) is considered to be a rigid obstacle (see Fig. 5). The acoustic velocity V a and the acoustic pressure P are solutions of the linearized Euler equations with velocity c a = 340 m/s, density ρ a = 1.29 kg/m 3 and adiabatic compressibility coefficient µ a = 1/(ρ a c 2 a ), in the unbounded domain Ω = R 3 \ {ω ∪ Ω r } which excludes the rim and the plate:

     ρ a ∂V a ∂t + ∇P = 0 µ a ∂P ∂t + Div V a = 0 in Ω ( 22
)
Viscothermal losses in the air are ignored in the acoustic model. The normal component of the acoustic velocity vanishes on the rim: The coupling between the 3D sound field in Ω and the vibrating soundboard in ω obeys to the condition of continuity of the velocity normal components:

V a • n r = 0 on Ω r (23)
V a • e z = ∂u p ∂t on ω (24)
where e z completes the referential (e x , e y ) introduced for the describing the motion of the soundboard in ω ⊂ R 2 (see Fig. 5). Finally, the soundboard force f is the pressure jump:

[P ] ω = P | ω --P | ω + (25)
RR n°8181

where ω + and ω -stand for the both sides the plate. Again, the vibroacoustic system satisfies the following energy decay:

dE p,a dt ≤ 0 with E p,a (t) = E p (t) + E a (t) (26)
where E p (t) was defined in Eq. [START_REF] Fletcher | Normal vibration frequencies of a stiff piano string[END_REF], and

E a (t) = Ω ρ a 2 |V a | 2 dx dy dz + Ω µ a 2 P 2 dx dy dz. (27)

Piano model

In the complete piano model, the soundboard is coupled to the hammer-strings system, according to the description made in Section 2.4, and radiates in free space (see Section 2.5). Consequently, the force f in the system of equations ( 15) becomes:

f = - i F b,i (t) χ ω + [P ] ω . ( 28
)
At the origin of time, the hammer has an initial velocity, while all the other parts of the system are at rest. The global coupled system satisfies the energy decay property:

dE h,s,p,a dt ≤ 0 with E h,s,p,a (t) = E h,s (t) + E p (t) + E a (t), (29)
where E h,s (t) is defined in Eq. ( 14), E p (t) is defined in Eq. ( 18), and E a (t) is defined in Eq. ( 27). In our model, it is assumed that there is no dissipation in the strings -soundboard and soundboard -air coupling terms. In other words, all dissipative terms are intrinsic to the constitutive elements (hammer, strings, soundboard). The acoustic dissipation is consecutive to the radiation in free space, which corresponds to energy loss for the "piano" system.

Numerical formulation

We now turn to the discretization of the global piano model described in the previous sections. We have to solve a complex system of coupled equations, where each subsystem has different spatial dimensions, inducing specific difficulties. The hammer-strings part is a 1D system governed by nonlinear equations. The soundboard is a 2D system with diagonal damping. The acoustic field is a 3D problem in an unbounded domain. The selection of the appropriate numerical methods are governed by the necessity of constructing an accurate and a priori stable scheme.

The nonlinear parts of the problem (hammer-strings interaction, string vibrations), the coupling between the subsystems and, more generally, the size of the problem in terms of computational burden, requires to guarantee long-term numerical stability. In the context of wave equations, and in musical acoustics particularly [START_REF] Rhaouti | Time-domain modeling and numerical simulation of a kettledrum[END_REF][START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF], a classical and efficient technique to achieve this goal is to design numerical schemes based on the formulation of a discrete energy which is either constant or decreasing with time. This discrete energy has to be consistent with the continuous energy of the physical system. If the positivity of the discrete energy is ensured, then a priori estimates can be established for the unknowns of the problem leading to the stability of the method. For most numerical schemes this imposes a restriction on the discretization parameters, expressed, for example, as an upper bound for the time step.

Another innovative aspect of our method is that the reciprocity and conservative nature of the coupling terms are guaranteed. In the discrete formulation, the couplings need a specific handling in order to guarantee a simple energy transfer without artificial introduction of dissipation, and without instabilities. Our choice here is to consider discrete coupling terms that cancel each other when computing the complete energy. In total, this method yields centered implicit couplings between the unknowns of the subsystems. The order of accuracy of the method is preserved, compared to the order of each subsystem taken independently, with no additional stability condition.

In view of the diversity of the various problems encountered in the full piano model, different discretization methods are chosen for each subsystem and for the coupling terms. The complete piano model is written in a variational form. In summary:

• Higher-order finite elements in space, and an innovative nonlinear three-points time discretization are used on the string,

• A centered nonlinear three time steps formulation is used for the hammer-strings coupling,

• A modal decomposition, followed by a semi-analytic time resolution is used for the soundboard,

• A centered formulation is used for the strings-soundboard coupling, where the forces acting at the bridge are introduced as additional unknowns,

• For the acoustic propagation, higher-order finite elements are used in the artificially truncated space, coupled to Perfectly Matched Layers (PML) at the boundaries, with an explicit time discretization.

The numerical schemes used are not presented in detail below. We restrict the presentation to some general survey on the numerical resolution and on its main difficulties. References are given to more numerically-oriented papers where a rigorous presentation of the method is given.

Strings

Standard higher-order finite elements are used for the space discretization of the nonlinear system of equations that govern the vibrations of the strings. The spatial discretization parameters (mesh size and polynomial order) are selected to ensure a small numerical dispersion in the audio range.

The unknowns are then evaluated on a regular time grid so that u(n ∆t) ≈ u n h , v(n ∆t) ≈ v n h and ϕ(n ∆t) ≈ ϕ n h . The most popular conservative schemes for wave equations are the family of θ-schemes, which have two drawbacks in the context of the piano strings: firstly, they were designed for linear equations, and, secondly, the less dispersive schemes in this family are subject to a restrictive stability condition that yields an upper bound for the time step (the so-called CFL condition). Therefore, it has been decided to adopt two different discretization schemes for the linear and the nonlinear part of the system, respectively. In addition, an improvement of the θ-scheme is developed for the linear part that combines both stability and accuracy.

For the nonlinear part, one major difficulty is that no existing standard scheme has the capability to preserve a discrete energy. As a consequence, we had to develop new schemes [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations. Application to the vibrating piano string[END_REF] based on the expression of a discrete gradient, which ensures the conservation of an energy, for a special class of equations called "Hamiltonian systems of wave equation". We have shown theoretically that no explicit scheme could ensure energy conservation, and the final numerical scheme is nonlinearly implicit in time (which implies that a nonlinear system must be solved at each time iteration). For a scalar equation, for example, the scheme would simply treat a continuous derivative term H ′ (u) as a derivative quotient based on the evaluation of the solution at successive time steps:

H ′ (u(n ∆t)) ≈ H(u n+1 h ) -H(u n-1 h ) u n+1 h -u n-1 h (30)
The case of a system of equations involves the gradient ∇H(u, v, ϕ) for which new methods have been developped since the trivial generalisation of Eq. ( 30) does not lead to an energy preserving discretisation.

In the linearized part of the string system, transverse, longitudinal and shear waves coexist. In the piano case, the transverse waves propagate much slower than the two others. In the audio range, and for the 88 notes of the instrument, small numerical dispersion must be guaranteed for the series of transverse partials and for the lowest longitudinal components. The shear components are beyond the audio range, and thus numerical schemes inducing higher dispersion for this series are acceptable.

In view of these considerations, a novel implicit discretization has been elaborated that reduces the numerical dispersion while allowing the use of a large time step in the numerical computations [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]. This method is based on the classical second-order time derivative centered scheme combined to the three points centered θ-approximation:

   ∂ 2 u ∂t 2 (n ∆t) ≈ u n+1 h -2u n h + u n-1 h ∆t u(n ∆t) ≈ θ u n+1 h + (1 -2θ) u n h + θ u n-1 h (31)
where θ ∈ (0, 1/2). For θ ≥ 1/4, the numerical scheme is unconditionally stable, while for θ < 1/4 the time step ∆t must satisfy the relation

∆t 2 ≤ 4 (1 -4θ) ρ (K h ) ( 32
)
where ρ (K h ) is the spectral radius of the stiffness matrix. When finite differences are used for space discretisation of a wave equation with velocity c, with a mesh size h, this relation simplifies to

c 2 ∆t 2 h 2 ≤ 1 1 -4θ . (33)
A classical analysis shows that the value θ = 1/12 minimizes the numerical dispersion, but choosing this value for the whole system would yield a too severe time step restriction, due to the two fastest waves. We propose a scheme where the value θ = 1/4 is used for the longitudinal and shear waves, hence relaxing the stability condition, while the value θ = 1/12 is used for the transverse wave, hence reducing the numerical dispersion for the series of transverse partials (see Fig. 6). The implicit nature of the resulting scheme might be a drawback in other contexts, but is not penalizing here since an implicit scheme is already necessary for the nonlinear part of the system. The solution is computed via an iterative modified Newton-Raphson method which needs the evaluation of both the scheme and its Jacobian with respect to the unknowns. It can be shown that a discrete energy is decaying, after extinction of the source. The stability of the numerical scheme

RR n°8181

Table 2: Parameters for the note D♯1 used in the simulations shown in Figure 6. The observation point is located at a distance of 6 cm from one end. Simply supported boundary conditions are considered. 2. Solid line: spectrum obtained from numerical simulation. Circles: theoretical spectrum of the numerical simulations. Diamonds: theoretical spectrum of the continuous system. (Top) Usual θ-scheme, with θ = 1/4. (Bottom) New scheme with θ = 1/4, θ = 1/12. is derived from this property, yielding a condition on the time step. In practice, for typical space discretization parameters, the time step ∆t = 10 -6 s yields stable and satisfying results in terms of dispersion.

L A ρ c F T 0 E (m) (m 2 ) (kg•m -3 ) (N) (Pa)

Hammer-strings coupling

The hammer-strings system is solved by considering together the unknowns of each strings of a given note, and the hammer displacement ξ n h ≈ ξ(n ∆t). The nonlinear hammer-strings interacting force is treated in a centered conservative way, following the method used for the strings system (the function Φ in Eq. ( 13) is seen as a derivative quotient of its primitive function -Ψ governed by Eq. ( 30)). A global discrete energy is shown to be decaying with respect to time when the hammer is given with an initial velocity.

Soundboard

The soundboard model assumes a diagonal damping in the modal basis. Its motion is first decomposed onto the modes of the undamped Reissner-Mindlin plate (see Eqs. ( 15)) belonging to the audio range, after semi-discretization in space with higher-order finite elements [START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF]. In practice, 2 400 modes are needed to model the soundboard vibrations up to 10 kHz. Fig. 7 shows some of the numerically computed modes associated to their eigenfrequencies. The presence of ribs and bridges is visible in the high frequency range. These modes are only computed once for all, before the time computation starts.

The choice of a Reissner-Mindlin model for the plate is validated numerically by comparison with measurements on rectangular plates made by other authors [START_REF] Lambourg | Modèle temporel pour la simulation numérique de plaques vibrantes. application à la synthèse sonore[END_REF] (see Fig. 8). The discrete modal amplitudes X h,n (t) of the discrete eigenmodes are solutions of the uncoupled second-order equations:

d 2 X h,n dt 2 + (2πf h,n ) 2 X h,n = F h,n
where the f h,n are the numerically computed eigenfrequencies of the undamped Reissner-Mindlin model. We then introduce a discrete damping mode per mode:

d 2 X h,n dt 2 + α(f h,n ) dX h,n dt + (2πf h,n ) 2 X h,n = F h,n (34)
This procedure yields decoupled equations which can be solved analytically in time, without introducing any additional approximation or numerical dispersion [START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF]. The energy identity over time of this semi-discrete problem is also exactly satisfied with this method. However, one drawback of this choice is the loss of the local nature of the coupling with strings and air.

Strings-soundboard coupling at the bridge

The discrete formulation of the strings-soundboard continuity equations must be done with special care in order to ensure the stability of the resulting coupled scheme. The purpose here is to couple the implicit three points nonlinear strings scheme described in Section 3.1 with the time semianalytic soundboard model described in Section 3.3. In practice, the last point of the string has to fulfil a discrete condition consistent with the continuity condition expressed in Eq. ( 19). For computational efficiency reasons, new variables are introduced that represent the coupling forces associated to the cinematic conditions between string and soundboard expressed in Eq. ( 19) (see Fig. 4). The strings and soundboard unknowns are evaluated on interleaved time grids: {n ∆t} for the strings, and {(n + 1/2) ∆t} for the soundboard. The forces at the bridge are considered constant on time intervals of the form [(n -1/2)∆t, (n + 1/2)∆t]. The coupling condition is an implicit version of Eq. ( 19) centered on times n ∆t (see Fig. 9).

Due to the linearity of the soundboard model, it is possible to express the soundboard unknowns as linear functions of the forces at the bridge. Taking advantage of this property, it is possible to perform Schur complements [START_REF] Brezinski | Schur complements and applications in numerical analysis[END_REF] on the system which, originally, is globally implicit. An algorithm is then written which updates first the unknowns of the strings and the forces at the bridge, and, in a second step, updates the unknowns of the soundboard.

Acoustic propagation and structural acoustics

The acoustic domain being unbounded, it is necessary to artificially truncate the computational domain, while minimizing the wave reflection on this artificial boundary. Perfectly Matched Layers is an option [START_REF] Bécache | Stability of perfectly matched layers, group velocities and anisotropic waves[END_REF]. Even after this truncation, the numerical parameters need to be chosen small enough, in view of the necessary wideband computation. A commonly accepted rule is to provide at least 10 points per wavelength so that the signal is spatially well sampled by the discretization. For example, if one has to model the propagation of an acoustic wave with noticeable energy content up to 10 kHz, with a corresponding smallest wavelength of 3.4 cm, one has to provide a mesh composed of 3 mm large elements. Since the size of a piano can be 2 m large and 3 m long, with a 40 cm rim, the mesh reaches 60 millions of degrees of freedom.

n n n + 1 n -1 n -1 n + 1 !"#$
h,n (t) X h,n (t) X h,n (t) P n+1/2 h P n-1/2 h V n a,h V n+1 a,h V n-1 a,h (u n h , v n h , ϕ n h ) (u n+1 h , v n+1 h , ϕ n+1 h ) (u n-1 h , v n-1 h , ϕ n-1 h )
Figure 9: Schematic view of the discretization. The string's variables (u h , v h , ϕ h ) are evaluated on the time grid n ∆t. The soundboard modal displacements X h,p are calculated at times (n+1/2) ∆t.

The acoustic velocity V a,h is calculated at times n ∆t and the acoustic pressure P h is calculated at times (n + 1/2) ∆t. All methods used yield energy identities, where the energies are centered on times (n + 1/2) ∆t. The coupling terms representing the forces at the bridge are centered on times n ∆t. (Color online)

The acoustical problem is solved in space with higher-order finite elements, and in time with an explicit leapfrog scheme, in view of the large number of degrees of freedom to consider. The acoustic velocity V a,h and sound pressure P h are calculated at times {n ∆t} and {(n + 1/2) ∆t}, respectively. This scheme has a restrictive stability condition: in practice, the adopted time step is around 10 -6 s. An implicit coupling exists between the soundboard displacement and the acoustic pressure in the vicinity of the plate, which implies a change of basis between both the physical and modal representations of the soundboard. In the variational formulation, the coupling between soundboard and air appears as skew-symmetric source terms for the soundboard and the sound pressure equations, respectively. These terms are constructed in the discrete scheme so that they vanish when computing the energy, and centered at times n ∆t. Due to the linearity of the equations, it is possible here to perform Schur complements, yielding an efficient algorithm that updates separately the plate (with a semi-analytic method) and the air variables.

A virtual piano

The resulting numerical scheme is stable, under the previously mentioned conditions on the numerical parameters, globally implicit, nonlocal for the soundboard, and uses multifarious methods. The efficiency of the computer code is optimized through the use of adapted additional unknowns and Schur complements, so that the update of the unknowns of each subsystem is made separately at each time step. A massively parallel computing was necessary, and special attention was paid to the cost of each step, in order to minimize the global computation time. In average, computing one second of sound for the complete piano model (with frequency content up to 10 kHz) takes 24 hours on a 300 cpus cluster (around 86 ms per time iteration for the complete piano).

In Fig. 10 the time evolution of selected quantities are represented for the note C2. Nearly 3 000 degrees of freedom (dofs) are used on each string of the triplet, 420 000 dofs are used to compute the 2400 first modes of the soundboard, and slightly more than 90 000 000 dofs are used in the air domain. The time step is 10 -6 second. The parameters are listed in Table 1 for the soundboard and in Table 3 for the strings. The hammer strikes the strings with a velocity of 4.5 m•s -1 (a forte to fortissimo playing). The longitudinal precursor can be seen on the upper parts of the figure: when the longitudinal wave of the string reaches the bridge, the soundboard is pushed down (see Fig. 10(b)), until the transverse wave arrives and pulls the soundboard up (see Fig. 10(g)). The acoustic wave is absorbed by the PML (not represented in the figures), and the rim is an obstacle to sound propagation (see Fig. 10(h)). Fig. 11 shows the energy evolution of each subsystem versus time, both in linear and logarithmic scales. The energy associated to sound propagation is computed in the truncated domain only. The total energy (solid line) is decreasing, as expected.

4

Results of the simulations

In order to show the ability of the method to simulate the complete register of a piano, the results (output data) of the model are presented and discussed in this Section for different notes in the low (D♯1 and C2), medium (F3) and treble (C♯5 and G6) register, respectively. Table 3 gives the numerical values used for these notes in the simulations. Since the nonlinear string is one major feature of the model, particular attention is paid to the influence of string amplitude on the resulting sound. The amplitude of transverse string motion is linked to the impact velocity of the hammer. The string amplitudes of the simulated tones, normalized with respect to string diameter, and the associated hammer impact velocities are shown in Table 4. Comparisons with real tones allows to group the simulations into typical categories of playing: piano, mezzo-forte, forte, fortissimo, corresponding to realistic, though relatively arbitrary, hammer velocity ranges [START_REF] Askenfelt | From touch to string vibrations: II: The motion of key and hammer[END_REF].

As a rule of thumb, it is generally admitted that, for thin structures, nonlinear effects become apparent when the amplitude of the vibrations are comparable to the thickness [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Transposing this rule to the case of strings yields the condition ADR = |u s | d s > 1, where d s is the diameter of the string. According to this rule, Table 4 predicts that nonlinear effects should be noticeable in the bass and medium ranges, even for moderate hammer velocity, whereas they should be less easily detectable for the treble notes, except for strong impacts. This simple criterion is in accordance with experimental observations. The order of the presentation follows the course of energy transmission in the piano: from hammer to strings, from bridge to soundboard, and from soundboard to air. The input parameter is the initial velocity of the hammer at the time where it comes in contact with the strings. The present model ignores the action of key mechanism prior to this contact. To illustrate this, Fig. 12 shows measured and simulated starting transients for the note D♯1. One can see the successive transformations of the hammer pulse to string wave, bridge and soundboard accelerations, and sound pressure. These waveforms and their corresponding spectra will be analyzed and discussed in the following.

An essential requirement for piano tone modeling is the accuracy of frequency estimation. To illustrate this feature, Fig. 13 shows an example of string inharmonicity for the note D♯1 (7th string with fundamental 39 Hz). Both measured and simulated string's eigenfrequencies follow the stiffness dispersion curve predicted by the Timoshenko model, at least up to the 60th partial (around 3 kHz). Precise measurements on partials of lower amplitude are difficult beyond this limit, because of noise and blurred spectral content.

One effect of string nonlinearity is the dependence of frequency with amplitude. For the note C♯5 played forte, for example, Fig. 14 shows that the frequency of the fundamental decreases with time, a consequence of amplitude decrease due to damping.

Observations made on real signals show that the pressure and soundboard motion spectra, including the bridge, have a denser and richer content than the strings. The simulations help here to understand these differences and identify the additional spectral components. In the low- frequency range, most of the additional spectral peaks correspond to soundboard modes excited by the string pulse. They are present in piano sounds even for light touch. These modes are particularly visible for the upper notes of the instrument, because of large spacing between the strings' partials (see Fig. 15). The soundboard modes are damped more rapidly than the string's partials, and thus they are essentially audible during the initial transients of the tones.

Increasing the hammer velocity progressively induces additional peaks between the string components, a consequence of string nonlinearity and coupling at the bridge. As explained in Section 2, a coupling exists between transverse and longitudinal motion of the nonlinear string. Due to string-bridge coupling, both the transverse and longitudinal components are transmitted to the soundboard. This explains why the longitudinal eigenfrequencies of the strings are visible on bridge and soundboard waveforms, but not on the transverse string motion [START_REF] Giordano | Motion of a piano string: Longitudinal vibrations and the role of the bridge[END_REF]. In addition, a consequence of the nonlinear terms in the string wave equation is that combinations of string components are created: the so-called "phantom partials" [START_REF] Conklin | Generation of partials due to nonlinear mixing in a stringed instrument[END_REF].

As explained theoretically from Eq. 21 in Section 2, the frequencies of these phantoms correspond to sums (or differences) between two or three of the string components, depending on whether the combinations are the results of quadratic, or cubic, nonlinearities. In general all combinations are not observable, and obey to complex rules. The instability conditions that gives rise to these components must reach a certain threshold [START_REF] Nayfeh | Nonlinear oscillations[END_REF]. Similar phenomena are observed in gongs and cymbals [START_REF] Touzé | Lyapunov exponents from experimental time series: Application to cymbal vibrations[END_REF]. Such an analysis is beyond the scope of this paper. Notice in Fig. 16 that the frequencies predicted by sums or differences of partials' eigenfrequencies correspond with great accuracy (less than 1 Hz) to the observed phantoms, both in simulations and measurements.

In the time-domain, the presence of nonlinear coupling is seen in the precursors. Fig. 12 shows that the precursor is not visible on the string displacement, while it is clearly seen on the other waveforms (bridge, soundboard, pressure), both in measurements and simulations. This precursor is primarily due to the longitudinal waves. In real instruments, the composition of these precursors is more complex and not totally understood [START_REF] Askenfelt | Observations on the transient components of the piano tone[END_REF]. Presumably, they contain some signature of shock of the key on the keybed, and other structural components that are not yet included in our model.

One further interest of simulations lies in the possibility of separating phenomena that are mixed together in the reality. In this respect, damping factors of the partials are good illustrating examples. When systems are coupled, it is always problematic to separate the causes of losses. In contrast, a model has the capability of introducing dissipation of energy in the hammer felt, along the string, and at the ends, separately. Even more interesting is the separation of structural and radiation losses. Experimentally, such a separation would require a rather delicate procedure where the instrument should to be put in a vacuum chamber in order to modify the conditions of radiation. One known difficulty of such experiments follows from the modifications of wood properties consecutive to variations of ambient pressure, and thus investigating such a problem with the help of simulations is an appealing alternative. To illustrate this ability of the model, Fig. 17 shows the damping factors derived from simulations for the C♯5 note, introducing each cause of losses successively (radiation, soundboard, strings, hammer felt). These damping factors are averaged in the frequency band of each partial, during the first second of the simulated tones. The damping factor related to the 17th partial (around 7 kHz) is ignored, since its amplitude is very low due to the striking position. In this example, the influence of the soundboard seems to be weak. For some partials, it turns also paradoxically out that the mean damping factors is less in the presence of both soundboard losses and radiation than for radiation only. In fact, it might be plausible that the conditions of coupling between strings and soundboard modes are slightly modified by the damping due to eigenfrequency shift. As a consequence, a "local" measure of damping might exhibit unexpected results. More investigation is needed here based, for example, on the computation of sound power and radiation efficiency.

Conclusion

In this paper, a global model of a grand piano has been presented. This model couples together the hammer, the nonlinear strings, the soundboard with ribs and bridges, and the radiation of acoustic waves in free field. As far as we are aware, this is probably the most general physical model of a piano available today. However, a number of significant features of real pianos were not considered in this model. The key mechanism, and an accurate description of the hammer action including the vibrations of the hammer shank, have been left aside. From the point of view of the player, adding this part would allow interesting studies of the links between the action of the finger and the resulting sound. Beside this, an improvement in the string model would be to account for the nonplanar motion observed on real pianos. A strong hypothesis is that this motion might be due to the customary observed zig-zag clamping conditions at the bridge [START_REF] Anderssen | The challenge of piano maker[END_REF], but this needs to be verified and quantified by measurements. To reproduce such effects, more appropriate boundary conditions have to be developed, that allow the progressive transformation of a vertical polarization into an horizontal motion. The way both polarizations are transmitted to the rest of the structure also need to be better understood. The coexistence of these polarizations with different decay times greatly influence the amplitude envelope of the tone, and thus its perception.

The motion of the structure is restricted here to the soundboard. Previous measurements tend to show that some other parts of the instrument contribute to the sound [START_REF] Askenfelt | Observations on the transient components of the piano tone[END_REF]. In this context, it would be attractive to reproduce the shock of the key against the keybed and the vibrations of the rim, to evaluate their relevance.

The present model is solved in the time-domain. The results yield the temporal evolution of the main significant variables of the system simultaneously: hammer force, string motion, bridge and soundboard vibrations, pressure field. The obtained waveforms can be heard through headphones or loudspeakers, and clearly evoke piano tones. They also shed useful light on the transfer of energy and transformations of the signals from hammer to strings, soundboard and air.

To numerically solve the problem, specific and original methods were developed for each part of the piano. A gradient approach coupled to higher-order finite elements lead us to design an energy decaying numerical scheme for the string's system, which is nonlinearly coupled to the hammer. A modal method was chosen to solve the soundboard problem, with a diagonal damping form. The eigenmodes and eigenfrequencies are computed once for all using higher-order finite elements, and an analytic formula is then used in time. Finally, sound radiation is solved with finite differences in time and higher-order finite elements in the space domain, which is artificially truncated with PML (Perfectly Matched Layers).

As a result, a numerical formulation of the global piano model is obtained with high precision in time, space and frequency. This formulation ensures that the total energy of the system is decaying. The model accounts for the dependence of piano sounds and vibrations with amplitude, due to nonlinear modeling of strings and hammers. In this respect, the simulations show the main effects of nonlinearity observed on real tones: precursors, time evolution of eigenfrequencies, transverselongitudinal coupling, and phantom partials. The model used for string-soundboard coupling at the bridge is consistent with the transmission of nonlinearities observed on real instruments. Due to this coupling, the presence of the soundboard modes in the piano transients are reproduced in a natural way. The soundboard model also integrates the presence of ribs and bridges, which are treated as heterogeneities in material and thickness of a Reissner-Mindlin plate.

The piano is a instrument with a large register. Most of the notes, from bass to treble, show a wideband spectrum, with significant energy up to 10 kHz and more. As a consequence, piano modeling requires a fine spatial grid for each part. The most demanding grid is associated with the modeling of the pressure field. This part of the simulations was a particularly challenging point of the work, which necessitated high performance parallel computing. In the present state of the equipment, several hours of computation in parallel on a 300 cpus cluster are necessary to compute the pressure field during one second in a box of the order of 10 m 3 that contains the instrument [START_REF] Chabassier | Physical parameters for piano modeling[END_REF].

Analysis of the simulated piano tones in time and frequency show a satisfactory agreement with measurements performed on a Steinway D grand piano. This particular instrument was used for extracting accurate values of input parameters, thus allowing precise comparisons between model and measurements for some selected notes in the bass, medium and treble range. Informal auditory evaluation of the simulated tones indicates that medium and treble notes are fairly well reproduced, but that the depth of the bass notes is not completely rendered.

In its present state, this model of piano should be considered as a crude skeleton of the instrument. Its prime function is to get a better understanding of the complex coupled phenomena involved in a complete piano, with the possibility of systematic variations of making parameters. Numerous additional improvements, careful adjustments and fine tuning would be necessary before thinking of competing with high-quality pianos. However, even in its imperfect form, we believe that the model could be used as a companion tool for piano making. In this context, investigating the influence of soundboard modifications on the radiation of sound and on string-bridge coupling appear as potentially fruitful examples.

Figure 1 :

 1 Figure 1: Schematic view of the string motion with the three unknowns: flexural displacement u s (x, t), longitudinal displacement v s (x, t), shear angle ϕ s (x, t).

Figure 2 :

 2 Figure 2: Spreading function δ H (x) used to model the hammer string contact. The width of the contact zone is about 2 cm.

Figure 3 :

 3 Figure 3: Schematic view of the four different zones of the soundboard in the referential plane. The orthotropy axis of the soundboard makes an angle of -40 degrees with the horizontal axis (see the stripes).

Figure 4 :

 4 Figure 4: (a) Schematic view of strings-soundboard coupling at the bridge. The soundboard is flat because of the static action of the strings. The bridge is supposed to move in the vertical direction ν only. The string forms a small angle α with the horizontal plane containing the vector τ . (b) The spot indicates the spreading function χ ω for note C2, centered on the point where the string passes over the bridge.

Figure 5 :

 5 Figure 5: Geometrical configuration of the piano. (Color online)

Figure 6 :

 6 Figure 6: Spectrum of the transverse displacement of string D♯1 when considering the linear Timoshenko model. Parameters are listed in Table 2. Solid line: spectrum obtained from numerical simulation. Circles: theoretical spectrum of the numerical simulations. Diamonds: theoretical spectrum of the continuous system. (Top) Usual θ-scheme, with θ = 1/4. (Bottom) New scheme with θ = 1/4, θ = 1/12.

Figure 7 :Figure 8 :

 78 Figure 7: Examples of soundboard eigenmodes, computed with fourth-order finite elements (nearly 450 000 degrees of freedom). Coordinates are in m. (Color online)

  %&%'()*+,-./0%1232$4'%& 50#$%&%'()*+,-*67#$%389 :&679'%&();<-.673=>62?=()@(<-.'?%38()A(<-9#/&'?2$(B%8B(6?=/?(+!C C6=2$(=/&60#69%'%63 D%8B(6?=/?(+!C

X

Figure 10 :Figure 11 :

 1011 Figure 10: 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 = 0.59 m and y = 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. (Color online)

Figure 12 :

 12 Figure12: Measured (left) and simulated (right) starting transients of the main variables for note D♯1 (7th). From top to bottom: hammer acceleration, string displacement (at point located 1.749 m from the agraffe), bridge acceleration at string end, soundboard acceleration (at point x=0.17 m ; y=1.49 m in the coordinate axes shown in Fig.3), sound pressure (simulated at point x=0.8500 ; y=1.4590 ; z=0.3800, and measured in the nearfield at a comparable location).

Figure 13 :

 13 Figure 13: Inharmonicity derived from frequency analysis of simulated (black circles) and measured (squares) string spectra. Note D♯1.

Figure 14 : 2 Figure 15 :Figure 16 : 1 )Figure 17 :

 1421516117 Figure 14: Evolution of frequency with time of the fundamental, due to geometrical nonlinearity. Simulation of note C♯5 played fortissimo.

Table 3 :

 3 Parameters used for the strings in the simulations -p ) 4.0 × 10 8 2 × 10 9 1.0 × 10 9 2.8 × 10 10 2.3 × 10 11 × 10 11 Pa, G = 8 × 10 10 Pa, ρ c = 7850 kg.m -3 , ρ = ρ c F.

		Note	D♯1	C2	F3	C♯5	G6
	L	(m)	1.945	1.600	0.961	0.326	0.124
	d	(mm)	1.48	0.9502	1, 0525	0, 921	0.794
	F		5.73	3.55	1	1	1
	T 0	(N)	1781	865	774	684	587
	f 0	(Hz)	38.9	65.4	174.6	555.6	1571.4
	p		2.4	2.27	2.4	2.6	3.0
	K H (N•m M H (g)	12.00	10.20	9.00	7.90	6.77
	x H	(m)	0.25	0.2	0.115	0.039	0.015
	x 0	(m)	0.47	0.59	0.54	0.88	1.16
	y 0	(m)	1.63	1.26	0.83	0.23	0.05
	A = , E = 2.0 RR n°8181 πd 2 4 , I = πd 4 64

Table 4 :

 4 String amplitude and initial hammer velocity of the tones simulated with the present model. ADR is the maximum amplitude-to-diameter ratio of the strings.

	String ADR Hammer velocity (m/s) dynamics
	D♯1	0.568	0.5	p
	D♯1	1.77	1.5	mf
	D♯1	3.59	3	f
	C2	0.685	0.5	p
	C2	2.01	1.5	mf
	C2	4.05	3	f
	F3	0.337	0.5	p
	F3	1.04	1.5	mf
	F3	2.10	3	f
	C♯5	0.274	0.5	p
	C♯5	0.87	1.5	mf
	C♯5	1.76	3	f
	C♯5	2.64	4.5	f f
	G6	0.12	0.5	p
	G6	0.443	1.5	mf
	G6	0.96	3	f

RR n°8181

Avenue de la Vieille Tour, 33405 Talence Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgment

Simulations 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/) and 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 (see http://www.mcia.univ-bordeaux.fr).