Juliette Chabassier

Antoine Chaigne

Patrick Joly

Time domain simulation of a piano. Part 1 : model description

Keywords: piano, modeling, energy, nonlinear precursor, damping phenomena * Magique 3d team piano, modélisation, énergie, précurseur non linéaire, phénomènes d'amortissement

The purpose of this study is the time domain modeling of a piano. We aim at explaining the vibratory and acoustical behavior of the piano, by taking into account the main elements that contribute to sound production. The soundboard is modeled as a bidimensional thick, orthotropic, heterogeneous, frequency dependant damped plate, using Reissner Mindlin equations. The vibroacoustics equations allow the soundboard to radiate into the surrounding air, in which we wish to compute the complete acoustical field around the perfectly rigid rim. The soundboard is also coupled to the strings at the bridge, where they form a slight angle from the horizontal plane. Each string is modeled by a one dimensional damped system of equations, taking into account not only the transversal waves excited by the hammer, but also the stiffness thanks to shear waves, as well as the longitudinal waves arising from geometric nonlinearities. The hammer is given an initial velocity that projects it towards a choir of strings, before being repelled. The interacting force is a nonlinear function of the hammer compression. The final piano model is a coupled system of partial differential equations, each of them exhibiting specific difficulties (nonlinear nature of the string system of equations, frequency dependant damping of the soundboard, great number of unknowns required for the acoustic propagation), in addition to couplings' inherent difficulties.

Introduction

This work is the continuation of a long term collaboration between the Unité de Mécanique, ENSTA ParisTech, specialized in musical acoustics and the project team POems (CNRS / ENSTA /INRIA) specialized in the development of numerical methods for wave equations. This collaboration, whose aim is to desing physical models and perform time domain simulations of musical instruments, already gave birth in the past to modeling tools for musical instruments, the timpani [START_REF] Rhaouti | Time-domain modeling and numerical simulation of a kettledrum[END_REF] and the guitar [START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF]. By considering today the piano, we attack a new challenge that represents a gap compared to the two above mentioned work as well as for the complexity of the physics and the underlying models as for the size of the problem.

This work is of course is related to the problem of sound synthesis, whose one aim is to generate realistic sounds of given instruments (here, the piano). Many methods reach this goal successfully. A number of them operate in real-time, based on various strategies (pre-recording of some selected representative sounds ; frequency, amplitude or phase modulation ; additive or subtractive synthesis ; parametric models...), see [Gui]. However, most of them have only little connections with the physics of the instrument.

Our wish is not only to reproduce the sound generated by a physical object (the piano) convincingly, but rather to understand how this specific object can generate such a particular sound, by modeling the complete instrument, based on the equations of the physics and using geometric and material related coefficients. Such an approach can be referred to as « physics based sound synthesis ».

The first works in this direction rely on very simplified or reduced models. With regard to the specific case of the piano, one can mention, for example, the very popular method of digital waveguides, through which a piano model was proposed by [START_REF] Nakamura | Piano tone synthesis using digital filters by computer simulation[END_REF] and reviewed by [BAB + 03]. This method is particularly effective and efficient and can be coupled with more traditional methods (as for instance finite differences in [START_REF] Bensa | The simulation of piano string vibration: from physical models to finite difference schemes and digital waveguides[END_REF]). An alternative approach, based on the so-called modal method, has been chosen for instance at IRCAM for the software Modalys [107].

Another approach is to use the standard tools of numerical analysis for PDEs (finite element, finite differences) to solve the system of equations numerically. The advantage of such an approach is to keep a strong connection to the physical reality, and to make very few a priori assumptions on the behavior of the solution. The intention is to reproduce the attack transients and the extinction of the tones faithfully, thus offering a better understanding of the complex mechanisms that take place in the vibrating structure. This approach was adopted in the past to study separate parts of the piano: [START_REF] Boutillon | Model for piano hammers: Experimental determination and digital simulation[END_REF] and [START_REF] Chaigne | Numerical simulation of piano strings. i. a physical model for a struck string usin finite-difference methods[END_REF] investigated the interaction between the hammer and the strings. [START_REF] Bilbao | Conservative numerical methods for nonlinear strings[END_REF] was interested in modeling the nonlinear behavior of the strings. [START_REF] Cuenca | Modélisation du couplage corde -chevalet -table d'harmonie dans le registre aigu du piano[END_REF] proposed a model to explain the coupling between the strings and the soundboard at the bridges. [START_REF] Vyasarayani | Modeling the dynamics of a compliant piano action mechanism impacting an elastic stiff string[END_REF] and [START_REF] Izadbakhsh | Dynamic modeling and experimental testing of a piano action mechanism with a flexible hammer shank[END_REF] studied the vibration the hammer shank. To our knowledge, there is only one published work [START_REF] Giordano | Physical modeling of the piano[END_REF] which focuses on both the modeling a full piano and its numerical formulation. This model makes use of partial differential equations, and accounts for the main involved phenomena : from the initial blow of the hammer to the propagation of sound, including the linear vibration of strings and soundboard. For the discrete formulation of the problem, classical numerical analysis tools are used, namely finite differences in space and time.

Our work aims at continuing this effort by providing a complete piano model (to our knowledge, it is the most accurate model available today), and a reliable, innovative and accurate numerical method to solve it. These two steps were naturally the two main parts of the work that has given rise to this series of two articles, the first of which is devoted to the construction of the mathematical model. The second part will be concerned by the discretization of this model and the validation of physical hypotheses through numerical simulations. Although this first article is supposed to be readable independently, its interest will appear more clearly when reading the second article. The goal of the present article is threefold:

(i) Explain the historical construction of our piano model, pointing out the links with the physics and the limitations of this model, (ii) Describe some fundamental mathematical properties (in particular energy identities) that provide a real confidence for its soundness from a theoretical point of view, (iii) Propose a general and abstract framework and formulation for this model, which we find useful for at least two reasons Inria -besides the fact that it permits some conciseness in the presentation, this prepares the second part on the discretization of the model, which will rely in an essential way on this formulation, -we believe that future enrichments of our model will fit this general framework, which should limit the amount of work for the improvement of our computational code.

Before giving the outline of this paper, it is first useful to describe in some detail the structure of a piano, introducing the main elements that we shall systematically refer to, as well as the main physical mechanisms involved in the sound production. These are illustrated by figures 1(a) and 1(b). The keyboard of most pianos has 88 black and white keys, corresponding to the notes of the tempered scale. The action of one key (see figure 1(a)) throws one hammer toward one or several strings, depending on the selected note. The strings are made of steel, but the bass strings are wrapped with copper. Each string is attached to a wooden beam, the bridge, which transmits its vibrations to the soundboard, a thin (less than 1 cm) wooden plate which radiates in the surrounding air inducing our perception of a sound. A cast iron frame is placed above the soundboard in order to support the strings' tension, and the complete system is integrated in a thick structure with the keyboard in front. Most of the time, three pedals are at the pianist's feet's disposal, allowing to act on the dampers or the hammers' mechanism. This principle is identical for grand and upright pianos, even though the practical implementation is different in each case. In the following, we adopt the Anglo-Saxon notation (from A to G) to name the different notes of the piano, starting from octave 0. The first key is therefore A0, corresponding to a fundamental frequency of 27.5 Hz, while the last one is C8, with fundamental 4186 Hz.

In the following of this paper, in section 1, we present various experimental results that will serve for assigning some objectives to our model. We especially describe some features that seem to be (at least all together) specific of a piano sound such as the sound precursor, the inharmonicity, and the phantom partials. In the long section 2, we present the model we have designed for the piano strings (a visco-elastic nonlinear stiff string model), for which a special effort has been given. We explain in particular why the simple linear model that we used for guitar strings in [START_REF] Derveaux | Time-domain simulation of a guitar: Model and method[END_REF] could not be satisfactory. Section 3 is concerned with the model we have chosen for the hammers and their interaction with the strings. Section 4 is devoted to the mathematical model for the soundboard (subsection 4.1) and more importantly for its coupling with the strings at the bridge (subsection 4.2). Finally, in section 5, we construct our full piano model by putting together the models on the previous sections and coupling them to the model for sound radiation described in section 5.1. Let us notice that a particular attention has been given to accounting for various damping mechanisms, which appear to be essential in sound perception.

Nowadays, the piano is certainly the most widely played instrument, and its advanced manufacture keeps up with its popularity. Beyond the aspects that concern fundamental research, which we have already emphasized, we hope that this work could also be of some interest for piano makers. The challenges they face today include the seek for volume and homogeneity from bass to treble, and even more : seek for a specific timbre (or tone color), long sustain, and for an appropriate distribution of sound in space. In order to reduce the proportion of empiricism, and anticipate the impact of possible changes in the vibrational and acoustical behavior of the instrument, many piano makers have their own research laboratory, oriented towards experimentation but also towards numerical simulation, as for instance the case of house Schimmel. Numerical methods are now part of the improvement and testing process of various parts of the piano (soundboard modal analysis, spectral analysis of the strings, shape optimization of the cast iron frame...). Some piano makers collaborate with universities or research laboratories, in order to answer specific questions about the instrument (radiation efficiency, characteristics damping time, or boundary conditions at the bridge as, for example, in the collaboration between the pianos Stuart & Sons and the Australian research centre CSIRO). The approach used by piano makers, however, suffers from one major limitation : although they are able to study in detail the behavior of each part of the instrument, they generally do not consider the coupling between its main elements, which in fact may significantly influence this behavior. In contrary, a comprehensive modeling tool, as we intend to design in this paper, accounting for all the couplings between the main parts of the instrument, yields a better understanding of the influence of some particular settings on the whole behavior of the piano. It becomes then possible to conduct « virtual experiments », by systematically changing materials, geometries, or some other design parameters, and observe the effect of these changes on the entire vibro-acoustic behavior of the instrument, and ultimately, on the resulting sound.

1 Experimental observations : sound precursors, inharmonicity and phantom partials

Experimental results and objectives for our model

We begin by a presentation of some experimental facts that are taken from the literature or are the results of measurements performed by two of the authors (J. Chabassier and A. Chaigne). These observations have served to define some objectives for the construction of our model.

We begin by analyzing a real piano sound which is nothing but a pressure signal p(t) as a function of time as illustrated in figure 2 perform a time frequency analysis by computing seismograms that result from windowed Fourier transforms of p(t)

p(t, f ) = +∞ 0 p(τ ) θ τ -t T e -2iπf t dτ (1)
Inria where θ(σ) is a normalized, smooth enough, window function satisfying

   θ ≥ 0, supp θ = [0, 1], θ > 0 in ]0, 1[ θ is increasing on [0, 1/2], θ(1 -σ) = θ(σ), θ(1/2) = 1 (2)
and T denotes the size of the window. The spectrogram is then defined as

S p(t, f ) = 20 ln p(t, f ) p ∞ , p ∞ = sup t,f | p(t, f )| (3)
In principle, T is chosen sufficiently large with respect to the periods of the oscillations of the signal and sufficiently small with respect to its duration. The spectrograms that we present in level lines in the (t, f ) plane (see figure 2(b)) have been computed with T = 200ms and the window function θ is the one of the MATLAB package. Roughly speaking, in each spectrogram :

• the curves f → S p(t, f ) shows the frequency content of the signal in [t, t + T ],

• the curves t → S p(t, f ) shows the time evolution of the component of the signal with frequency f .

The two spectrograms of figure 2(b) correspond to the same note G3. The left one corresponds to a piano dynamic level play, the right one to a fortissimo dynamic level play. On each spectrogram, the presence of horizontal rays illustrates the existence of predominant frequencies and the variable decay of these rays illustrates that the attenuation of the signal is greater for high frequencies than for smaller ones.

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

If the two pressure signals were to each other, the two spectrograms would be identical, which is not the case. This indicates that some nonlinear phenomena have been involved in the physical mechanisms.

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

The next results that we wish to point out are extracted from the article [START_REF] Podlesak | Dispersion of waves in piano strings[END_REF] in JASA. The authors performed some experiments related to the note A0 (the first A of the keyboard) of a grand piano. They were interested in

• the transverse displacement (more precisely the "vertical" transverse displacement, parallel to the hammer striking direction) of a point of the string located 10 cm from the bridge,

• the sound pressure signal at a point corresponding to a microphone placed about 10 cm above the piano.

From a naive reasoning, one expects that the pressure signal begins after the string vibration. However, this is exactly the contrary that the authors observed as illustrated on figure 3(a), which concerns the first 80 ms of their experiment. The pressure signal begins with a "high frequency" (compared to what is observed on the transverse displacement) signal that clearly contributes to the sound : they called it the sound precursor.

They analyzed the frequency content of this part of the signal through its spectrum (in short, the modulus of the Fourier transform of the signal adequately truncated in time, represented in logarithmic scale) : this gives a curve with respect to the frequency f , where 2πf is the dual variable for the time, see figure 3(b) for the range of frequencies [0, 5000] Hz. One clearly sees regularly spaced peaks. These coincide with the ones observe on the longitudinal displacement of the string when is it excited only in its longitudinal direction.This in opposition with the usual transverse solicitation due to the hammer, which explains that such peaks cannot be observed on the spectrum of the transverse displacement. This was in some sense a "demonstration" that the longitudinal displacement has a role to play in the generation of a piano sound.

Objective 3. Account for some mechanism of transmission of the longitudinal string's displacements to the soundboard. From [START_REF] Podlesak | Dispersion of waves in piano strings[END_REF].

Next, we present some results of experiments that we did on a Steinway grand piano (D model) that was put at our disposal for some period at IRCAM1 . These experiments concern the note D 1. In figure 4(a), we observe, in the range of frequencies [1 2] kHz, the spectrum of the transverse displacement of the string at a point located 1.8 cm from the agraffe. Again we observe sharp peaks, indicated by the red circles, that correspond to an increasing sequence f n of frequencies. At first glance, these peaks seem to be equally spaced but a closer analysis indicates that the spacing between two consecutive peaks, f n+1 -f n , slightly increases with n. More precisely, these peaks seem to follow (at least in the range of audible frequencies [0, 20] kHz) a law of the form:

f n = n f 0 (1 + B n 2 ). ( 4
)
This is called the inharmonicity of a piano sound (as opposed to the case where peak frequencies would be all proportional to a fundamental frequency f n = n f 0 ). The dimensionless parameter B, which is small with respect to 1, is called the inharmonicity factor.

Objective 4. Reproduce the inharmonicity effects.

Finally, we observe on figure 4(b), the spectrum of the vertical2 acceleration of the bridge at the attached point of the string. We see, in addition to the peaks observed on figure 4(a), some new peak frequencies indicated by the magenta diamonds. These frequencies also appear on the spectrum of the vertical acceleration of any other point of the soundboard, as illustrated by figure 4(c), as well as in the spectrum of the recorded sound pressure signal, see figure 4(d). These partials were named "phantom partials" by Conklin in [START_REF] Conklin | Piano strings and "phantom" partials[END_REF] when he first observed them experimentally, since there was no explanation of there existence from existing models.

Objective 5. Account for the phantom partials.

A mathematical model for the piano string

For pedagogical purpose and for the sake of completeness, we are going to construct progressively our string model by successive modifications of the simplest possible model, namely the linearized vibrating string equation. We shall justify our successive enrichments of this model by the objectives of section 1. In our explanations, we do not pretend to a complete mathematical rigor (which would probably be out of reach or would demand lengthy and tedious technical developments) but aim at providing some intuition to the reader. When we will speak of linear or linearized models, we shall rather systematically refer to the spectral analysis of a "generalized harmonic oscillator". By this expression, we mean an abstract linear evolution equation of the form

d 2 U dt 2 + A U = 0, (5)
where the unknown function t → U (t) takes its values in some Hilbert space H and A denotes an unbounded positive selfadjoint operator in H, with compact resolvent. The notion of mode of vibration, or eigenmode, for equation ( 5) is linked to the research of a particular solution of the form

U (t) = U a e ±2iπf t , U a ∈ H, f ∈ R + (6)
which leads to find the frequency f such that λ = 4 π 2 f 2 is an eigenvalue of A. Denoting 4 π 2 f 2 n , n ≥ 0 the spectrum of A, where f n is an increasing sequence of real numbers that tends to +∞, the numbers f n are, by definition, the eigenfrequencies of the harmonic oscillator. Moreover, introducing an orthonormal basis W n ∈ H, n ≥ 0 of related eigenfunctions of A, it is well known that any "finite energy" solution of the free evolution problem (5) is of the form

U (t) = +∞ n=0 ± u ± n W n e ±2iπfnt (7)
where the real coefficients u ± n satisfy appropriate summability conditions. From the above formula, it is straightforward to establish the coincidence between the frequencies f n and the peaks of the spectrum of the "signal" U (t) (or a U T (t) e -2iπf t dt.

(8)

This will help us establish a link between a model of the form (5) and the experimental results of section 1.

u(x) x u s (x) (a) Transversal displacement x ϕ(x) ϕ s (x) u s (x) (b) Shear angle u(x) v(x) x u s (x) v s (x) (c) Longitudinal displacement

The vibrating string equation

The most common model for describing the transverse vibrations of a string assumes that the point of the string moves only transversely in a vertical plane along a line which is orthogonal to the string at rest (the reference configuration, represented by a segment [0, L]). As a consequence, the only unknown of the model in u s (x, t) the transverse displacement of the point of abscissa x (at rest) at time t (see figure 5(a)). Under the usual "small displacement" (u s remains small with respect to L) and "small deformations" (x-derivatives of u s remain much smaller than 1) assumptions, the displacement field is governed by the d'Alembert's equation (or 1D wave equation):

ρA ∂ 2 t u s -T 0 ∂ 2 x u s = 0, ∀ x ∈ [0, L], t > 0, (9)
where, assuming that the string is homogeneous and of constant cross-section, ρ denotes its volumic mass, A the area of its section and T 0 its tension at rest. The corresponding propagation velocity is

c τ = T 0 /ρA . ( 10
)
where the index τ refers to transverse.

If one completes (9) by homogeneous Dirichlet boundary conditions

u s (0, t) = u s (L, t) = 0, t > 0, (11)
which express that the string is fixed at its two extremities, one obtains, as it is well-known, an harmonic oscillator whose eigenfrequencies and corresponding eigenmodes w n are given by

f Ale n = nf 0 , f Ale 0 = c τ 2L , n ≥ 1, w n (x) = 2 L sin nπx L . ( 12
)
The set f Ale n , n ≥ 1 forms what is called a "harmonic spectrum", because it is made of frequencies that are the multiple of a fundamental frequency f Ale 0 which is the musical pitch of the note.

Remark 1. On a piano string, the Dirichlet condition at the first extremity of the string (let us say x = 0) is completely justified by experimental data. It is no longer true, of course, that the bridge extremity (x = L) does not move since this is precisely where the string's oscillations are transmitted to the soundboard. However, due to the strong rigidity of the bridge, the amplitude of these oscillations remains quite small in comparison with the ones of the string's central point. This is why the Dirichlet boundary condition can be seen as a reasonable approximate boundary condition for the string : in first approximation, analyzing the vibrations of the string in terms of the Dirichlet problem's eigenmodes provides reliable insights about the physical reality. We shall see how the condition (11) must in fact be modified when we shall treat the coupling with the soundboard (cf section 4.2).

Inria

Stiffness and inharmonicity : the linear stiff string equation

The prestressed Timoshenko's beam model.

Since the d'Alembert's equation is unable to represent the inharmonicity of a piano sound, a more elaborate model has to be used : the stiff string model. In this model, one does not only take into account the transverse displacement of the string but also the fact that each cross section of the string can rotate with respect to the normal plane to the so-called neutral fiber of the string, as illustrated by figure 5(b). That is why one introduces an angle ϕ s as a new unknown. One has then to take into account that a portion of string applies on the adjacent portions not only a tension but also a torque. The equations governing the variations of (u s , ϕ s ), called linear stiff string equations in the context of this paper, are given by the prestressed Timoshenko's beam model. These equations involve new geometric and material properties of the string, namely the inertia momentum of the string's cross section I and the Young's and shear moduli, namely E and G, of the material composing the string. The parameter κ, called "Timoshenko's parameter", or "shear correction factor" is a dimensionless parameter, between 0 and 1 (see [START_REF] Cowper | The shear coefficient in timoshenko's beam theory(shear coefficient in timoshenko beam theory)[END_REF] for a physical approach and [START_REF] Babuska | Optimal shear correction factors in hierarchical plate modelling[END_REF] for a mathematical discussion of its value).

     ρA ∂ 2 t u s -T 0 ∂ 2 x u s + AGκ ∂ x ϕ s -∂ x u s = 0, x ∈ ]0, L[, t > 0, (13a) ρI ∂ 2 t ϕ s -EI ∂ 2 x ϕ s + AGκ ϕ s -∂ x u s = 0, x ∈ ]0, L[, t > 0. (13b)
It is possible, via a "dimensional analysis" to interpret the Timoshenko model as a perturbation of the d'Alembert model, considering the small dimensionless parameter:

η := d L
where d is the diameter of the cross section of the string, (

where the string (at rest) is assimilated to a cylinder. Of course, the use of a 1D dimensional model is justified by the fact that η is small. From a dimensional analysis, we can write

A = η 2 A * , I = η 4 I * . (15)
On the other hand, to maintain the propagation velocity c τ of the transverse waves constant, and thus keep the same fundamental frequency f Ale 0 , one must compensate, in the limit process, the decay of the area of the cross section by decreasing the tension of the string, i.e. considering that T 0 = η 2 T * 0 , which implies that c τ and f Ale 0 (cf. ( 10) and ( 12)) are independent of η.

From the second equation of (13), one can write formally

∂ x ϕ s = 1 + η 2 (A * Gκ) -1 ρ I * ∂ 2 t -EI * ∂ 2 x -1 ∂ 2 x u s
which we can substitute into the first equation to obtain

ρA * ∂ 2 t u s -T * 0 ∂ 2 x u s + A * Gκ 1 + η 2 (A * Gκ) -1 ρ I * ∂ 2 t -EI * ∂ 2 x -1 -1 ∂ 2 x u s = 0 (17)
that leads to a fourth order (in space and time) partial differential equation for u after applying the operator ρ I * ∂ 2 t -EI * ∂ 2

x to both sides of the equation. However, the last term in ( 17) is small since formally:

1 + η 2 (A * Gκ) -1 ρ I * ∂ 2 t -EI * ∂ 2 x -1 -1 ∼ -η 2 (A * Gκ) -1 ρ I * ∂ 2 t -EI * ∂ 2 x (η → 0)
an this is why (17) is a ( second order in η) perturbation of the d'Alembert's equation.

In first approximation (in the sense of remark 1), equations (13) are naturally completed by "simply supported" boundary conditions, namely the homogeneous Dirichlet boundary conditions (11) for u s plus a condition of "zero torque" at each extremity of the string, which amounts to imposing homogeneous Neumann boundary conditions for

ϕ s : ∂ x ϕ s (0, t) = 0, ∂ x ϕ s (L, t) = 0, t > 0. ( 18
)
It is straightforward that (9, 11, 18) corresponds again to a generalized harmonic oscillator whose eigenfrequencies and corresponding eigenmodes can be determined analytically, thanks to the choice of boundary conditions (11, 18). We shall not detail here these analytical computations which are long and tedious but straightforward (see [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]) and shall restrict ourselves to describe the most useful results. As we have a system of two second order equations, it is not surprising that these modes can be split into two families of modes (the following splitting appears naturally in the analytical computations)

• the family of "flexural" modes, with frequencies f n , n ≥ 1 and eigenmodes (w n , ψ n ), n ≥ 1 ,

• the family of "shear" modes, with frequencies f S n , n ≥ 1 and eigenmodes

(w S n , ψ S n ), n ≥ 1 , that satisfy f n+1 > f n , f S n+1 > f S n , f S n > f n , ∀ n ≥ 1.
We represent in figure 6(a) the "curves" of both families of frequencies with physical data that correspond to the string A first striking fact is that shear modes eigenfrequencies are much larger than the flexural modes eigenfrequencies and are all above 20 kHz. As a consequence such modes can not contribute to a perceptible sound. Moreover, very high frequency sources should be needed for exciting these modes, which is not the case of the piano's hammer solicitations.

In figure 6(b), one gives a closer look of the first 200 flexural eigenfrequencies that all belong to the interval [0, 7500] Hz. These are represented by the blue diamonds that progressively deviate, when n increases from the black circles corresponding to the harmonic spectrum of the d'Alembert's equation. We observe that these frequencies are exhibit a behavior that is similar to the one of measured eigenfrequencies. To be more precise, using a Taylor expansion of their analytical expression gives:

f n ≈ nf 0 1 + B n 2 , B = π 2 EI T 0 L 2 1 - T 0 EA , for small enough values of n √ B ( 19
)
which correspond to the red triangles in figure 6(b). This shows that inharmonicity is indeed foreseen in this model, with an inharmonicity factor B (as defined in section 1, formula 4) which is typically, for real piano strings, of the order of 10 -4 . From the spectral point of view, the link between the Timoshenko and D'Alembert's model can be understood through an asymptotic analysis with respect to the small parameter η defined by (14). More precisely, one can show (modulo long and tedious computations) that, when η → 0

lim η→0 f n = f Ale n , lim η→0 f S n = +∞ in O(η -1 ) (20)
This gives a more analytic insight about what is observed on figures 6(a) and 6(b).

Remark 2. The fact that, in figure 6(b), the "curve" of the frequencies f n 's is above the "curve" of the frequencies f Ale n 's is due to the fact that, for real piano strings T 0 /EA << 1.

Timoshenko versus Euler-Bernoulli : a short discussion

There exists another well-known model that accounts for inharmonicity effects : the Euler-Bernoulli model. Contrary to the Timoshenko's model, the (scalar) unknown is the same than for the d'Alembert's model, namely the transverse Inria displacement u s . As the d'Alembert's model, the Euler-Bernoulli model can be recovered from a perturbation analysis for small values of the parameter η (see ( 14)), of the Timoshenko model. If we start from the equation (17) for u s , instead of dropping the last term as in the previous section, one can retain the following O(η 4 ) formal approximation

1 + η 2 (A * Gκ) -1 ρ I * ∂ 2 t -EI * ∂ 2 x -1 -1 ∼ -η 2 (A * Gκ) -1 ρ I * ∂ 2 t -EI * ∂ 2 x (η → 0).
Doing so, one obtains an approximate equation for u which is of second order in time:

ρ A * -η 2 I * ∂ 2 x ∂ 2 t u s -T * 0 -η 2 EI * ∂ 2 x ∂ 2 x u s = 0 (21)
One can get rid of the ∂ 2 x ∂ 2 t term by making the additional approximation

A * -η 2 I * ∂ 2 x ∼ A * , (22)
which leads to the so-called Euler-Bernoulli model which, going back to our original notation, takes the form

ρA ∂ 2 t u s -T 0 ∂ 2 x u s + EI ∂ 4 x u s = 0. ( 23
)
When one replaces the Timoshenko model by the approximate Euler-Bernoulli's model, one replaces a fourth order in time equation for u s (17) by a second order in time equation ( 23). As a consequence, the very high frequency modes associated to frequencies f S n disappear, leaving only one family of eigenmodes (as for the d'Alembert's equation) that are "approximations" (in η) of the flexural modes associated to frequencies f n of the Timoshenko's model. Nevertheless, contrary to d'Alembert's model, Euler-Bernoulli's model allows us to take into account inharmonicity effects since the corresponding frequencies are given by

f EB n = f Ale n 1 + n 2 π 2 EI T 0 L 2 1 2 (24)
which yields in particular

f EB n ≈ nf 0 1 + B EB n 2 , B EB = π 2 EI 2T 0 L 2 , for small enough values of n (25)
The reader will note that the Timoshenko and Euler-Bernoulli's model lead to different inharmonicity factors (cf. ( 19) and ( 25)). In practice, for real piano strings, since T 0 /EA is very small with respect to 1 (as already mentioned in remark 2), these are not much different. It is also difficult to decide which model is closer to the physical reality, even though the Timoshenko's model seems to be richer from the physical point of view. This question is still a subject of debate in the acoustics/mechanics community. The reason why we have preferred the Timoshenko's model in this work comes more from practical arguments : even it introduces an additional unknown, it avoids to deal with fourth order spatial differential operators, which is easier from the numerical point of view.

2.3 Taking into account the longitudinal displacement: the geometrically exact model

Derivation of the geometrically exact model

Up to now, we did not mention the longitudinal displacements of the strings. Let us forget for a while the rotations of the cross-sections of the string and the stiffness unknown ϕ s to come back to a model where the motion of the string is described only through the displacement of the neutral fiber's points, as for the vibrating string equation. The difference is that this displacement is allowed to have a longitudinal component, denoted v s (see figure 5(c)). This leads to the so-called "geometrically exact model" (GEM), as derived in [START_REF] Morse | Theoretical Acoustics[END_REF], that is established without referring to any "small displacement" or "small deformation" assumption. Let us denote u s (x, t) = u s (x, t), v s (x, t) the displacement vector at point x (in the reference configuration) and time t, and T(x, t) the tension of the string at the same point and same instant, which represents the action of the portion of string [x, L] on the portion [0, x]. From the fundamental law of dynamics, these are related by

∂ 2 t u s -∂ x T = 0. ( 26
)
The relation between T and the deformation of the string comes from the fact that • T(x, t) is tangent to the deformed string at (x, t) : T is collinear to (∂ x u, 1 + ∂ x v), which yields

T = T (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 ∂ x u, 1 + ∂ x v s t , T ∈ R, (27)
• T is the sum of the initial tension T 0 at rest and, by Hooke's law (assuming a linear behavior of the material), a term which is proportional to the infinitesimal elongation a(x, t) ∈ R of the string at point x and time t (meaning that a small element of string of length ∆x centered at point x at t = 0 has a length, after deformation, equal to (1 + a)∆x + O(∆x 2 ) ). This gives, thanks to elementary geometry

T = T 0 + E A a, a = (∂ x u s ) 2 + (1 + ∂ x v s ) 2 ) 1 2 -1. ( 28
)
Substituting the expression of T into (26) leads to the following 2 × 2 system of nonlinear equations :

                 ρA ∂ 2 t u s -∂ x EA ∂ x u s -(EA -T 0 ) ∂ x u s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 = 0, x ∈ ]0, L[, t > 0, (29a) ρA ∂ 2 t v s -∂ x EA ∂ x v s -(EA -T 0 ) (1 + ∂ x v s ) (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 = 0, x ∈ ]0, L[, t > 0. ( 29b
)
Note that, since the material is assumed to have a linear behavior, the only nonlinearity of the model comes from geometrical nonlinearities (due to the elongation a) which justifies the name "geometrically exact" model.

In the case of a finite string, equation (33) has to be completed by boundary conditions, for instance by expressing that the string is fixed at its two extremities (which is realistic in first approximation for the piano, see remark 1 again) which gives (11) for u s and

v s (0, t) = v s (L, t) = 0, t > 0, (30)

Mathematical structure and properties of the model

One can check that this system can be put into an second order hamiltonian form by introducing the elastic potential energy density function H(u, v) : R 2 → R + defined by

H(u, v) = 1 2 EA (u 2 + v 2 ) -(EA -T 0 ) u 2 + (1 + v) 2 1 2 -(1 + v) (31)
Verifying that H is a positive function is left to the reader. Is is then easy to verify that the tension T of the string is given by

T = ∇H(∂ x u s , ∂ x v s ) (32)
which means that (29) can be rewritten as:

ρA ∂ 2 t u s -∂ x ∇H(∂ x u s ) = 0. (33)
One can then show that (33) also enters the category of (locally) nonlinear hyperbolic systems. Indeed, introducing the new unknown vector:

U = (∂ t u s , ∂ x u s ) ∈ R 4 (34)
(33) can be rewritten as a first order system

∂ t U + ∂ x F (U) = 0. ( 35
)
where, denoting U = (U t , U x ) ∈ R 2 × R 2 the current vector of R 4 , the flux function F is given by

F (U ) ≡ F (U t , U x ) = -(ρA) -1 ∇H(U x ), U t t . ( 36
)
The Jacobian of F has the following 2 × 2 block decomposition Figure 7: Level sets around the point (0, 0) of the energy density H(u, v) and its approximations. The physical parameters have been chosen so that T 0 /EA = 0.6, which emphasizes the visualization of the convexity loss. fact that DF (U ) is diagonalizable with real eigenvalues at least for |U | small enough, follows from the local (strict) convexity of H(u, v). This is deduced from the Taylor expansion of H(u, v) around the origin

DF (U ) = -(ρA) -1   0 D 2 H(U x ) I 0   (
H(u, v) = H 2 (u, v) + O (u 2 + v 2 ) 2 , H 2 (u, v) = T 0 2 u 2 + EA 2 v 2 . ( 38
)
We can visualize the region of convexity of H in the (u, v) plane (and thus the region of hyperbolicity of (35)) on figure 7(a) where we represent the level lines of H(u, v). In addition, it can be shown that the system is, in its region of hyperbolicity, linearly degenerate. That is, if ± λ j (U ), j ∈ {1, 2} are the (real) eigenvalues of DF (U ) with corresponding eigenvectors (in R 4 ) r ± 1 (U ), r ± 2 (U ) , then ∇λ j (U ) ∈ R 4 is orthogonal to r ± j (U ), j ∈ {1, 2} (the calculations are done in [START_REF] Chabassier | Energy preserving schemes for nonlinear hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF]). This has nice mathematical consequences, which seem to be physically relevant : in particular, the existence and uniqueness of smooth (C 2 in space and time) global solutions for the Cauchy problem associated to (33) (or equivalently (35)) provided that the initial data are smooth enough, in the C 2 norm (see for instance [START_REF] Ta-Tsien | Global classical solutions for quasilinear hyperbolic systems[END_REF]).

An important consequence of the structure of (33) is that sufficiently smooth solutions satisfy an energy conservation result, namely

d dt E s (u s ) = 0, E s (u s ) = 1 2 L 0 ρA ∂ t u s 2 dx + L 0 H(∂ x u s ) dx. ( 39
)
This provides a fundamental stability property for the model since a priori estimates (in H 1 -norm) are easily deduced from this energy conservation result.

The linearized model

In the case of small amplitude motions, it is natural to use a linearized model obtained by replacing in (33) H(u, v) by its quadratic approximation H 2 (u, v) (see (38)). Doing so, one obtains two decoupled 1D wave equations, the first of which coincides with the vibrating string equation (9):

     ρA ∂ 2 t u s -T 0 ∂ 2 x u s = 0, with velocity c τ = T 0 /ρA (40a) ρA ∂ 2 t v s -E A ∂ 2 x v s = 0, with velocity c = E/ρ (40b)
where the index stands for longitudinal. Taking realistic values for ρ, A, E and I for real piano strings one observes that c c τ ≥ 10,

which means that longitudinal waves propagate much faster than transverse waves : this explains the role of the string's longitudinal vibrations in the existence of the sound precursor. However, for our purpose, the decoupled linear model (40) is not satisfying. Indeed in the case of a transverse solicitation as the hammer's, a source term will appear only at the right hand side of the first equation, which means that v s will remain identically 0 and that the longitudinal vibrations cannot be observed. This objection no longer holds for the exact model (29) because of the nonlinear coupling between the two equations. Even if a source term appears only on the first equation of , v s will

Chabassier & Chaigne & Joly not remain 0 since T (u, 0) = 0 ! This is the first motivation to keep the nonlinear model.

Nevertheless, the linear model (40) will be useful to "analyze" (in first approximation) the behavior of the solution of the exact model (29) in the case where the deformations are not too large : it gives the tangent harmonic oscillator at the origin of the nonlinear evolution problem (29). In particular, the spectrum of this harmonic oscillator is made of the union of two harmonic spectra:

• the transverse harmonic spectrum

f Ale n = n f Ale 0 , n ≥ 1 , f Ale 0 = c τ 2L
• the longitudinal harmonic spectrum

f n = n f 0 , n ≥ 1 , f 0 = c 2L
Note that in practice, because of the velocity contrast (41), f n ≥ 10 f Ale n and the transverse spectrum is much denser than the longitudinal spectrum. However the intersection between longitudinal spectrum and the range of audible frequencies, namely, [0, 20] kHz ∩ f n = n f 0 , n ≥ 1 is mostly non empty for all piano strings : for instance, for the note D 1, in this interval of frequencies, there are typically 20 "longitudinal eigenfrequencies" (versus 200 "transverse frequencies" ).

An approximate model with polynomial nonlinearities

Still in the case of small amplitude motion, one can obtain approximate models (hopefully more accurate than the linear (40)) by replacing H(u, v) by other approximations at the origin that H 2 (u, v). Let us present below a model, which presents the a priori advantage to generate only polynomial nonlinearities, which is one of the reason why, for instance, it has been used model in [START_REF] Bilbao | Conservative numerical methods for nonlinear strings[END_REF] for numerical approximation (see also [START_REF] Bank | Generation of longitudinal vibrations in piano strings: From physics to sound synthesis[END_REF] for more analytical purposes). This model is obtained from an "anisotropic" quartic approximation of H(u, v) in the sense that it is obtained from a Taylor expansion which is of fourth order in u but only second order in v:

H(u, v) = H app (u, v) + O(u 4 + v 2 ), H app (u, v) = 1 2 T 0 u 2 + T 0 2 v 2 + EA -T 0 2 v + u 2 2 2 . ( 42
)
This type of anisotropic approximation is justified in the case where the string is transversally solicited (see [START_REF] Chabassier | Transitoires de piano et non linéarités des cordes : mesures et simulations[END_REF]).

To be more explicit let us introduce the two functions

T τ (u, v) = ∂H ∂u (u, v), T (u, v) = ∂H ∂v (u, v) (43)
so that T τ (∂ x u s , ∂ x v s ) and T (∂ x u s , ∂ x v s ) are respectively the transverse and longitudinal components of the tension T of the string and consider the equations with a transverse source term, of small amplitude ε, hence for the equation in u s only:

   ρA ∂ 2 t u ε s -∂ x [T τ (∂ x u ε s , ∂ x v ε s )] = ε f (t), (44a)
ρA ∂ 2 t v ε s -∂ x [T (∂ x u ε s , ∂ x v ε s )] = 0, (44b)
It it easy to see formally that

u ε s = O(ε), v ε s = O(ε 2 ). (45)
In other words longitudinal vibrations have a much smaller amplitude than transverse ones. As a consequence, from (42), one deduces that

H(u ε s , v ε s ) = H app (u ε s , v ε s ) + O(ε 4 ). ( 46
)
According to (42), one has

         T τ (u, v) = T 0 u + (EA -T 0 ) uv + 1 2 u 3 ) + O u 4 + v 2 (47a) T (u, v) = EA v + 1 2 (EA -T 0 ) u 2 + O u 4 + v 2 (47b)
so that, replacing H by H app in (33), we obtain the following coupled system of equations

           ρA ∂ 2 t u s -∂ x T 0 ∂ x u s + (EA -T 0 ) ∂ x u s ∂ x v s + 1 2 (∂ x u s ) 3 = 0, (48a) ρA ∂ 2 t v s -∂ x EA ∂ x v s + EA -T 0 2 (∂ x u s ) 2 = 0, (48b)

Inria

In our work, we shall use the exact model but this model will be helpful for the interpretation of some of the numerical results provided by our model.

Remark 3. Proceeding as in section 2.3.2, it is clear that the first order system corresponding to the approximate model ( 48) is still locally hyperbolic (see also figure 7(c)). However, it can be shown that it is genuinely nonlinear. As a consequence "shocks", namely discontinuities of ∂ t u s and ∂ x u s (which implies the presence of kinks in the deformed shape of the string, which seems unphysical) will develop in finite time, even with arbitrarily smooth and small data.

The energy conserved for smooth solutions of ( 48) is of course

E app s (u s ) = 1 2 L 0 ρA ∂ t u s 2 dx + L 0 H app (∂ x u s ) dx. ( 49
)
The reader will note that the positivity of H app (u, v), and thus good stability properties of the model via energy estimates, is only granted provided that EA ≥ T 0 , which was not needed for the exact model but is nevertheless true for real piano strings.

Combining longitudinal vibrations and inharmonicity: the nonlinear stiff string model

The model for planar motions

The model we shall propose for modeling a piano string aims (in the absence of damping phenomena -see section 2.5) at combining the inharmonicity effects obtained with the linear stiff string model of section 2.2 with the longitudinal / transverse vibrations coupling effects provided by the geometrically exact model of section 2.3. That is why we proposed a model with three scalar unknowns:

• the transverse component of the displacement (in a vertical plane) : u s (x, t),

• the longitudinal component of the displacement : v s (x, t),

• the angle of rotation of cross sections (in a vertical plane) : ϕ s (x, t), that is obtained by "concatenation" of both models (13) and (29). Seen as a modification of the geometrically exact model, this consists in the following operations:

• In the equation (26) for the displacement field u s = (u s , v s ), the tension T, given by ( 27) for the geometrically exact model (29), is modified by adding the contribution due to the rotation of the cross sections, namely (note that only the transverse component of the tension is affected)

T = T (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 (∂ x u, 1 + ∂ x v s ) t + AGκ (ϕ s -∂ x u s ), 0 t ( 50
)
where T is still given by (28).

• the two equations for (u s , v s ) are completed by the equation governing ϕ s from the Timoshenko model.

As a consequence, one computes that the longitudinal and transverse component of the tension T are given by

             T τ = EA ∂ x u s -(EA -T 0 ) ∂ x u s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 + AGκ (ϕ s -∂ x u s ) T = EA ∂ x v s -(EA -T 0 ) (1 + ∂ x v s ) (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 (51) RR n°8097
Chabassier & Chaigne & Joly which leads to the following nonlinear system that we shall refer to as the nonlinear stiff string model:

                               ρA ∂ 2 t u s -∂ x EA ∂ x u s -(EA -T 0 ) ∂ x u s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 + ∂ x AGκ (ϕ s -∂ x u s ) = 0, x ∈ ]0, L[, t > 0, ρA ∂ 2 t v s -∂ x EA ∂ x v s -(EA -T 0 ) (1 + ∂ x v s ) (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 = 0, x ∈ ]0, L[, t > 0, ρI ∂ 2 t ϕ s -∂ x EI ∂ x ϕ s + AGκ ϕ s -∂ x u s = 0, x ∈ ]0, L[, t > 0.
(52)

This system will be completed, in first approximation, by the boundary conditions (11, 30, 18). Note that the "tangent" harmonic oscillator to this system at the origin is made of the (decoupled) union of the Timoshenko's model ( 13) in (u s , ϕ s ) (with boundary conditions (11 -18) ) with the 1D wave equation (40b) for v s (with boundary conditions (30)

). The spectrum of the linearized model is thus the union of three parts

f n , n ≥ 1 ∩ f n , n ≥ 1 ∩ f S n , n ≥ 1 (53)
where we recall that (section 2.2)

• the first part f n , n ≥ 1 is an inharmonic spectrum corresponding to flexural modes (section 2.2),

• the second part f n , n ≥ 1 is a harmonic spectrum corresponding to longitudinal modes (section 2.3),

• the third part f S n , n ≥ 1 (shear modes) does not intersect in practice the set of audible frequencies.

We shall call the model (52)the stiff nonlinear string model. For conciseness of our presentation, it will be useful to rewrite (52) in a more compact and abstract form. This is the object of the next subsection.

An abstract model for generalized non linear string equations.

Let us start from an energy density function of 2N variables:

H(p, q) : R N × R N -→ R + . ( 54
)
We shall denote respectively ∇ p H(p, q) ∈ R N and ∇ q H(p, q) ∈ R N the partial gradients of H(p, q) with respect to p and q respectively. Moreover, splitting the set of indices {1, • • • , M } as:

{1, • • • , M } = I D ∪ I N , I D ∩ I N = ∅, (55)
we introduce the orthogonal projector J from R N onto

V D = q = (q j ) 1≤j≤N | q j = 0 if j ∈ I D . (56)
The abstract string model then reads

                           Find q(x, t) : [0, L] × R + -→ R N such that M ∂ 2 t q -∂ x (∇ p H(∂ x q, q)) + ∇ q H(∂ x q, q) = 0 x ∈ ]0, L[, t > 0, J q(0, t) = J q(L, t) = 0, t > 0.
Id -J ∇ p H(∂ x q, q)(0, t) = Id -J ∇ p H(∂ x q, q)(L, t) = 0, t > 0.

(57) in which q(x, t) ∈ R N is the vector of "generalized" string unknowns while ∇ p H(∂ x q, q) is the vector of generalized efforts from which we can define the generalized tension

T = J ∇ p H(∂ x q, q) (58)
Inria which coincides with the physical tension T given by (50) when H is given by (63)) while

(I -J) ∇ p H(∂ x q, q)
is the "generalized" torque. The boundary conditions can be interpreted as "mixed Dirichlet-Neumann" boundary conditions in the sense that :

• {q j , j ∈ I D } are the unknowns to which a homogeneous Dirichlet condition is applied,

• The conditions ∂ pj H(∂ x q, q) = 0, j ∈ I N , are generalized homogeneous Neumann conditions.

Once again, an energy conservation result is associated with (57). More precisely, any sufficiently smooth solution of (57) satisfies:

d dt E s (q) = 0, E s (q) = 1 2 L 0 M ∂ t q • ∂ t q dx + L 0 H(∂ x q, q) dx. ( 59
)
The proof is quite standard and the details are left to the reader. One takes the inner product in R N of the first equation of (57) and integrate in space over [0, L]. Then, the following two ingredients are used for treating the second term of the first equation of (57), after integration by parts,

• For the integral part, one uses the chain rule

∂ t H(∂ x q, q) = ∇ p H(∂ x q, q) • ∂ 2 xt q + ∇ q H(∂ x q, q) • ∂ t q. (60)
• For the boundary terms at x = 0 and L, one writes

∇ p H(∂ x q, q) • ∂ t q = (Id -J)∇ p H(∂ x q, q) • ∂ t q + ∇ p H(∂ x q, q) • J ∂ t q. (61)
It is an exercise to check that (52, 11, 18, 30) enters this general framework with:

N = 3, I D = {1, 2}, q = (u s , v s , ϕ s ), M =   ρA 0 0 0 ρA 0 0 0 ρI   (62)
and the energy density function given by

H(p, q) = 1 2 T 0 |p 1 | 2 + 1 2 EA |p 2 | 2 + 1 2 EI |p 3 | 2 + 1 2 AGκ |q 3 -p 1 | 2 + (EA -T 0 ) 1 2 |p 1 | 2 + (1 + p 2 ) -p 2 1 + (1 + p 2 ) 2 . ( 63
)
For numerical purposes again, it will be useful to separate the energy density function H(p, q) as the sum of its "quadratic part", defined as the quadratic form which approaches H(p, q) at third order at the neighborhood of the origin, from its "non quadratic part", namely the rest. Looking more closely at the expression of H in (63), we see that it is of the form (note in particular that the non quadratic part of H only depends on p):

H(p, q) = H 2 (p, q) + U (p), U (p) = O(|p| 3 ) (64)
where H 2 (p, q) is given by ( • represents below the inner product in R N or R 2N depending on the context) With ( 64) and ( 65), the partial differential equation in (57) can be rewritten

H 2 (p, q) = 1 2   A B B t C     p q   •   p q   ≡ 1 2 A p • p + C q • q + 2 B q • p (
M ∂ 2 t q -∂ x A ∂ x q -∂ x (B q) + t B ∂ x q + C q + ∂ x (∇U (∂ x q)) = 0 (66) Chabassier & Chaigne & Joly
In (66), ∂ x (∇U (∂ x q)) represents the "nonlinear part" of the model while the other terms constitute the "linear part" of the model or equivalently the "tangent" harmonic oscillator. These two parts of the model will be treated differently when we shall deal with the time discretization of the problem.

In the particular case of the energy density (63), (64, 65) hold with

A =   T 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κ   , (67)
U (p) = (EA -T 0 ) p 2 1 2 + (1 + p 2 ) -p 2 1 + (1 + p 2 ) 2 . ( 68
)
Note that, separately, H 2 (p, q) and U (p) are not necessarily positive (for instance, this is not the case with (67, 68)) but their sum is. Such a property is in particular satisfied if one can decompose the matrix A as

A = A S + A T such that C t B B A
T is a positive matrix, and

1 2 A S q • q + U (q) ≥ 0, ∀ q ∈ R N . ( 69
)
In the following, we shall assume that (69) holds. This will be important for the stability analysis of our numerical method. In particular, (69) holds for (67, 68)) with:

A S =   T 0 0 0 0 EA 0 0 0 EI   . ( 70
)

An enriched string model authorizing non planar motions

It is possible to enrich the string's model (52) while remaining in the general framework (57). It is in particular possible to take into account the so-called double polarization of the string, which amounts to authorizing non planar motions. This leads to introduce a second transverse component for the string displacement, orthogonal to the preponderant one u s . Proceeding in the same way for the rotation of the cross sections of the string, it is natural to introduce two additional unknowns (whose meaning is given in figure 8):

• u s : the second (or horizontal) component of the transverse displacement of the string,

• ϕ s : the second (or horizontal) angle for the rotations of the cross-sections, Doing so, we obtain a model (we do not write the equations in detail, these are straightforward extensions of (52) ) of the form (57) with :

N = 5, I D = {1, 2, 3}, q = (u s , u s , v s , ϕ s , ϕ s ).
Various experimental studies show that piano strings do have horizontal movements under the "vertical" solicitation of the hammer and that taking into account the double polarization of the string may have some importance from the acoustical point of view. However, we shall not consider them in the rest of this paper. This would be relevant only if we worked with a hammer model explaining how to generate horizontal vibrations of the string and with a bridge model explaining how these vibrations are transmitted to the soundboard. This will not be the case of the "simplified" models that we shall consider (see sections 3 and 4) and this is why we shall restrict ourselves to the "planar model" of section 2.4.1.

The full nonlinear stiff string model with intrinsic damping

It seems essential to account for the frequency dependent damping on strings. Damping phenomena are difficult to apprehend, for many reasons (lack of measurements, uneasy dissociation of origins, misunderstanding of certain phenomena as dislocation. . . ). This is why we propose as a first approach to use a very simple model that allows to feign this effect without trying to model the underlying physics. We have chosen to mimic the introduction of damping terms in d'Alembert's equation ( 9) by simply adding viscoelastic terms: where R u and γ u are empirical (constant in space and positive) damping coefficients which are respectively homogeneous to a time or the inverse of a time.

ρA ∂ 2 t u s + 2ρA R u ∂u s ∂t -2 T 0 γ u ∂ 3 u s ∂t ∂x 2 -T 0 ∂ 2 x u s = 0, ∀ x ∈ [0, L], t > 0, ( 71
) Inria u s (x) v s (x) ũs (x)
By extension of (71), we have chosen to treat our string's model ( 52) in a similar way by adding linear viscoelastic damping terms to each line of the system:

                                             ρA ∂ 2 t u s + 2 ρA R u ∂ t u s -2 T 0 γ u ∂ 2 x ∂ t u s -∂ x EA ∂ x u s -(EA -T 0 ) ∂ x u s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 + AGκ ∂ x ϕ s -∂ x u s = 0, ρA ∂ 2 t v s + 2 ρA R v ∂ t v s -2 T 0 γ v ∂ 2 x ∂ t v s -∂ x EA ∂ x v s -(EA -T 0 ) 1 + ∂ x v s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 = 0, ρI ∂ 2 t ϕ s + 2 ρI R ϕ ∂ t ϕ s -2 EI γ ϕ ∂ 2 x ∂ t ϕ s -∂ x EI ∂ x ϕ s + AGκ ϕ s -∂ x u s = 0, (72)
where (R u , R v , R ϕ ) and (γ u , γ v , γ ϕ ) are heuristic positive damping coefficients whose value is determined in practice thanks to experimental calibration. Again, in first approximation, this system is completed by the boundary conditions (11, 30, 18).

For conciseness and sake of generality (see section 2.4.3), we shall put the above model in an abstract and concise framework, using the notation of section 2.4.2 (see (54), ( 57)):

M ∂ 2 t q + ∂ t R q -∂ x (Γ ∂ x q) -∂ x (∇ p H(∂ x q, q)) + ∇ q H(∂ x q, q) = 0 ( 73
)
where R and Γ are N × N positive and symmetric matrices representing the damping terms. In the particular case where H(∂ x q, q) is of the form (64) -( 65), this gives

M ∂ 2 t q + ∂ t R q -∂ x (Γ ∂ x q) -∂ x (∇ p H(∂ x q, q)) + ∇ q H(∂ x q, q) = 0 (74)
According to section 2.4.2, the reader will easily check that (72) corresponds to the abstract form (74) where N = 3, q = (u s , v s , ϕ s ), the matrices M, A, B and C are given by ( 62) and (67), the function U is given by (68), while the

R = 2   ρA R u 0 0 0 ρA R v 0 0 0 ρI R ϕ   , Γ = 2   T 0 γ u 0 0 0 EA γ v 0 0 0 EI γ ϕ   , (75)
The reader will easily check that the boundary conditions (11, 30, 18) correspond to the following ones for the abstract model (where I D = {1, 2}, see ( 55) and ( 56))

J q(x, t) = 0, Id -J Γ ∂ 2 xt q + ∇ p H(∂ x q, q) (x, t) = 0, x = 0 or L, t > 0. ( 76
)
Moreover, it is important to notice that, in the presence of damping, the (generalized) tension of the string is no longer given by (58) but by

T = J ∇ p H(∂ x q, q) + Γ ∂ 2 xt q. ( 77
)
which gives in particular, if (62) holds, H is given by ( 63) and Γ by (75) and q by , a 2d vector whose longitudinal and transverse components are (instead of (51))

             T τ = EA ∂ x u s -(EA -T 0 ) ∂ x u s (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 + AGκ (ϕ s -∂ x u s ) + γ u T 0 ∂ 2 xt u s , T = EA ∂ x v s -(EA -T 0 ) (1 + ∂ x v s ) (∂ x u s ) 2 + (1 + ∂ x v s ) 2 1 2 + γ v T 0 ∂ 2 xt v s (78)
In this case, it is immediate to check that the energy conservation (59) for ( 57) is replaced by an energy decay result (that emphasizes the roles of R and Γ as damping coefficients)

d dt E s (q) + L 0 R q • q dx + L 0 Γ ∂ x q • ∂ x q dx = 0, E s (q) = 1 2 L 0 M q • q dx + L 0 H(∂ x q, q) dx, (79)
where we denote ẋ the time derivative of any variable x.

3 A mathematical model the strings-hammers interaction. At first glance, a piano's hammer (represented by H in what follows) can be described as a non deformable piece of wood covered by a deformable piece of felt. Each hammer will interact with one or several strings: for most notes, strings are gathered into "choirs" of one, two or three parallel strings that contribute to the same note. In what follows Inria we shall denote N s the number of strings. For each string, we shall use the model with planar motion of section 2.4.2. We will denote the i th string's unknowns q i = (u i , v i , ϕ i ), 1 ≤ i ≤ N s (for simplicity of the notation, we omit the index s for the string's unknowns) and by x ∈ [0, L] the abscissa along each of these strings. The strings of a same choir are slightly detuned (their tension at rest, T 0 , is different, see [START_REF] Weinreich | Coupled piano strings[END_REF]) and this is why, thanks to (67, 68), each string has its own U i and A i in (64, 65). Describing in detail the physics of the interaction between the hammer and the related strings (as a 3D contact problem for instance) would lead to a too complex model. For now, we shall restrict ourselves to a (very) simplified model.

The first part of the model is the kinematic one. In first approximation, the movement of the hammer is supposed to be parallel to a line D, vertical and orthogonal to the plane containing the strings with which it interacts. For simplicity:

• the movement of the wooden part of the hammer will be described by the abscissa ξ(t) (along D) of a fixed point M H in this wooden part (see figure 9(a)),

• the deformation of the felt will be described by the abscissa along D, ξ i (t), of the point of impact between the hammer and the i th string of the choir (which is also supposed to move along a line parallel to D, that we assume to be oriented in such a way that ξ(t) < ξ i (t))

• the above impact point is supposed to coincide with a point of abscissa x i along the string (we suppose that it does not depend on time : there is no slipping) so that

ξ i (t) = u i (x i , t) (80)
If one wants to take into account that the contact string-hammer is not purely punctual but distributed along a small portion of string around the point x i (doing so, we implicitly assume that the zone of impact does not change in time, which is also a simplification)

ξ i (t) = L 0 u i (x, t) δ H (x -x i ) dx (81)
where δ H (x) is a function with small support around the origin and satisfying δ H (x) = 1 .

To be more precise, if ξ > 0 denotes the distance between M H and the upper par of the hammer when this one is at rest (see figure 9(a))

there is contact with the i th string ⇐⇒ ξ i (t) -ξ(t) < ξ (82)
Moreover we shall define the crushing (see figure 9(b)) of the hammer at point x i by

e i (t) = ξ + ξ i (t) -ξ(t) + (83)
where x + is the positive part of x. In particular e i (t) = 0 when there is no contact.

The second part of the model consists in describing the interacting forces between the hammer and the strings. More precisely, we assume that the action of the hammer on the i th string is reduced to a force, applied at the impact point x i (or more specifically distributed along a small portion of string around x i , see later). We suppose that this force is orthogonal to the string and thus described by a scalar function F i (t). When one neglects the relaxation phenomena (see below), the intensity of this force is supposed to be a function of the crushing

F i (t) = k H Φ H e i (t) (84)
where k H > 0 is the nonlinear stiffness of the hammer (see [START_REF] Stulov | Hysteretic model of the grand piano hammer felt[END_REF]) and Φ H : R → R + is increasing and satisfies Φ H (0) = 0, so that, by construction

• when there is no contact (e i (t) = 0), this force vanishes,

• when there is contact, the more compressed the hammer is, the more intense the force is.

The reader will remark that the discrimination between k H and Φ H will only make sense when we shall have chosen an explicit expression for Φ H , see (86). For a more realistic modeling, one should take into account the relaxation phenomenon that expresses that the force of interaction is not the same depending on the fact that the hammer is being compressed, in which case the force is more intense, or decompressed, in which case it is less intense. Mathematically, this can be translated into the addition to (84) of a term which is proportional to the time derivative of Φ H e i (t) :

F i (t) = k H Φ H e i (t) + r H d dt Φ H e i (t) (85)
where r H ≥ 0 is the relaxation coefficient. The model (85) allows to account for the hysteretic behavior of the hammer observed in experimental studies as [START_REF] Stulov | Hysteretic model of the grand piano hammer felt[END_REF] and is also responsible, as we shall see, for dissipation phenomena. The same studies also show that it is sufficient to take a simple form for the function Φ H :

Φ H (y) = y p (86)
where the exponent p, that depends on the considered hammer, is not necessarily an integer and in practice varies between 1.5 and 3.5.

In conclusion, introducing the vector ν s = (1, 0, 0) t as an element of the q space R N , so that u i = q i • ν s , the evolution of the hammer-strings system will be modeled by the following set of equations

                           Find ξ(t) : R + -→ R and q i (t) : [0, L] × R + -→ R N , 1 ≤ i ≤ N s such that m H d 2 ξ dt 2 (t) = - i∈Is 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 (x, t) • ν s δ H (x -x i ) dx, 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 (87)
where the last equation holds for each 1 ≤ i ≤ N s . The above system is completed by the boundary conditions (76) for each q i and by the initial conditions:

       ξ(0) = -ξ, dξ dt (0) = v H , q i (x, 0) = ∂q i ∂t (x, 0) = 0, x ∈ [0, L], (88)
which express that the string is initially at rest and that, at time t = 0, the hammer strikes the strings from below with an initial velocity v H . In (87), we shall clarify that

• Φ H is given by (86), m H > 0 is the effective mass of the hammer, k H and r H its stiffness and relaxation coefficients, the functions H i (p, q) are defined as in section 2.4.2 (through (63) or equivalently (63, 67, 68)) and the the matrices (M, R, Γ) as in section 2.5, except for the apparition of the index i, devoted to express the dependence with respect to the string, of the various coefficients of the string model, especially the initial tension.

• the presence of the function δ H , in the right hand side of the last equation, allows to represent the spatial distribution of the force exerted by the hammer. For energy conservation reasons, this has to be the same function as the one in the definition of ξ i (t). The presence of the vector ν s = (1, 0, 0) t (as an element of R N , the space of q unknowns) means that the direction of the force is along the line D, which implies that the right hand side only impacts the equation for u i .

The above model is associated to an energy decay result. More precisely, the energy of the strings-hammer system is defined by

E s,h (q i , ξ) = i E s (q i ) + m H 2 ξ 2 + i k H Ψ H e i , e i (t) = ξ + ξ i (t) -ξ(t) + , ξ i (t) = ξ i (t) = L 0 q i (x, t) • ν s δ H (x -x i ) dx (89)
Inria where the function Ψ H is defined by

Ψ H (x) = x 0 Φ H (y) dy (≥ 0). (90)
In (89), in addition to the strings' energies (the first term), one identifies the kinetic energy of the hammer (the second term) and the energy of interaction between the strings and the hammer (the last term). It is easy to establish the following energy identity satisfied by any (smooth enough) solution of (87),

d dt E s,h (q i , ξ) = - i L 0 R qi • qi - i L 0 Γ ∂ x qi • ∂ x qi - i r H Φ H e i ∂ t e i 2 (91)
Note that, since Φ H is increasing, (91) illustrates the role of the relaxation term in (85) as a dissipation term.

4 A mathematical model for the soundboard -strings interaction.

4.1 The mathematical model for the soundboard

The Reissner-Mindlin model

The soundboard's thickness being very small compared to the other two dimensions, we have chosen to model this structure as a plate (the index p in what follows will refer to "plate") on a bidimensional domain ω. We shall denote x = (x 1 , x 2 ) ∈ R 2 the space variable on ω (and z the third "vertical coordinate" according to figure 10). Seen as a 3D object, ω is in the plane z = 0. In the piano context, it is reasonable to consider only small displacements and small deformations. We use a standard linearized model, the Reissner-Mindlin model [START_REF] Reissner | The effect of transverse shear deformation on the bending of elastic plates[END_REF] which can be seen as a 2D version of the Timoshenko model ( 13) for stiff strings. The unknowns of the model are (see also figure 10) • the transverse displacement u p of the plate,

• the two deflection angles (θ 1,p , θ 2,p ) = θ p representing the rotations of the normal fibers of the plate.

The coefficients appearing in the model are

• the thickness δ > 0 the plate, its volumic mass ρ > 0 ,

• the elasticity coefficients : the shear and stiffness tensor G and C (2 × 2 positive definite symmetric matrices, diagonalizable in the same basis)

• a shear correction coefficient of the model denoted κ 2 > 0

RR n°8097

The tensor C permits to take into account the orthotropy of the wood. The thickness will be taken as a function of x, which allows to take into account the diaphragmatic nature of the soundboard (the fact that the thickness of the plate varies slowly and smoothly, from 6 cm in the middle to 9 cm at the boundaries, [START_REF] Conklin | Design and tone in the mechanoacoustic piano. part ii. piano structure[END_REF]) but also the presence of the bridge and the ribs, seen on figure 11. The other coefficients will also be varying with x which accounts for any possible heterogeneity of the materials. If the plate is submitted to a surfacic density f (x, t) of transverse forces, the governing equations are the following:

         ρ(x) δ(x) ∂ 2 t u p -div δ(x) κ(x) 2 G(x) ∇u p + θ p ) = f (x, t), ρ(x) δ(x) 3 12 ∂ 2 t θ p -Div δ 3 (x) 12 C(x) ε(θ p ) + δ(x) κ 2 (x) G(x) ∇u p + θ p = 0 in ω (92)
where we have used the following notation

• Div σ denotes the (vectorial) divergence of a 2D tensor field σ:

Div σ i = 2 j=1 ∂ xj σ ij , i = 1, 2
• ∇θ (resp. ε(θ)) is the tensor corresponding to the gradient (resp. its symmetric) of a 2D vector field θ :

∇θ ij = ∂ xj θ j , i, j = 1, 2, ε(θ) = 1 2 ∇θ + ∇θ t
Of course, (92) has to be completed by boundary conditions. Roughly speaking, the plate ω is "fixed" along its boundary ∂ω. However, the optimal boundary conditions are very difficult to estimate in realistic configurations. We have considered ideal simply supported conditions, namely (n denotes the unit normal vector to ∂ω): Remark 4. By anticipation, let us make precise that, in the full piano model, the transverse force f appearing at the right hand side of (92) will be the sum of two contributions f = f s + f a :

u p (x) = 0, ∇θ p (x) n = 0, ∀ x ∈ ∂ ω. ( 93
)
• f s : the force applied by the strings at the bridge (see section 4.2),

• f a ; the force applied by the outside air (see section 5).

Inria

The soundboard model in condensed form

In order to rewrite (92), (93), it is useful to introduce the vector of plate unknowns

U p = (u p , θ p ) t ∈ R M , M = 3, (94)
the diagonal matrix

M p (∇ x ) =       ρ(x) δ(x) 0 0 0 ρ(x) δ(x) 3 /12 0 0 0 ρ(x) δ(x) 3 /12       (95)
the differential operator

A p (x, ∇ x )U p = A p (x, ∇ x )   u p θ p   =    -div δ(x) κ(x) 2 G(x) ∇u p + θ p ) -Div δ 3 (x) 12 C(x) ε(θ p ) + δ(x) κ 2 (x) G(x) ∇u p + θ p    (96)
and the trace operator

B p (n, ∇ x )U p = B p (n, ∇ x )   u p θ p   =   u p ∇θ p n   . (97)
Then, denoting ν p = (1, 0, 0) t ∈ R M (the space for the U p unknown), (92) becomes

   M p (∇ x ) ∂ 2 t U p + A p (x, ∇ x )U p = f ν p , x ∈ ω, t > 0, B p (n, ∇ x )U p = 0, x ∈ ∂ω, t > 0. (98)
Again, there is a natural energy identity satisfied by any regular solution of (92), (93). More precisely, if we introduce the energy of the soundboard as

E p (U p ) = 1 2 ω ρ δ |∂ t u p | 2 + 1 2 ω ρ δ 3 12 |∂ t θ p | 2 + 1 2 ω δ 3 12 C ε(θ p ) : ε(θ p ) + 1 2 ω δ κ 2 • G |∇u p + θ p | 2 (99)
then, one has the identity

d dt E p (U p ) = ω f ∂U p ∂t • ν p (≡ ω f ∂u p ∂t ). (100)

Introducing intrinsic plate damping

As it was also the case for the string, there are intrinsic mechanisms, of various natures, that provoke damping phenomena, that is to say natural attenuation in time of plate vibrations. Describing in detail the physics of these mechanisms is a very complicated task that goes much beyond the objective of this work. This is why, as for the string, we have chosen to use a more heuristic and phenomenological model. A full thesis [START_REF] Ege | La table d'harmonie du piano -études modales en basses et moyennes fréquences[END_REF] has been recently devoted to the physics of soundboards and we shall use the results therein to design our damping model. In two short sentences, the conclusions of [START_REF] Ege | La table d'harmonie du piano -études modales en basses et moyennes fréquences[END_REF] (i) The soundboard damping is a linear phenomenon of modal nature, (ii) The attenuation increases with the frequency of the modes.

Mathematically, it is useful to reconsider the model (98) in the form of an harmonic oscillator (5) in the sense of the beginning of section 2 ( see (5)):

d 2 U p dt 2 + A p U p = 0, (101)
The Hilbert space H is H = H p = L 2 (ω) M equipped with the scalar product:

U p , U p H = M p (x) U p (x) • U p (x) dx (102)
Chabassier & Chaigne & Joly and A p U p = A p (x, ∇ x )U p for any U p ∈ D(A p ) where the domain of A p is

D(A p ) = {U p ∈ H 1 (ω) M / A p (x, ∇ x )U p ∈ L 2 (ω) M , B p (n, ∇ x )U p = 0 on ∂ω} (103)
The eigenmodes {W n : ω → R M , n ∈ N} and corresponding eigenvalues {λ n ≥ 0, n ∈ N} (which are ranked by increasing values and tend to +∞) of this harmonic oscillator are the solutions of the eigenvalue problem

   A p (x, ∇ x )W n = λ n M p (x) W n , x ∈ ω, B p (n, ∇ x )W n = 0, x ∈ ∂ω. (104)
After appropriate normalization, {W n , n ∈ N} form an orthonormal basis of the Hilbert space H. In order to introduce some dissipation in the model, while taking account of the observations (i) and (ii), we modify (101) in the following way:

d 2 U p dt 2 + f d A p dU p dt + A p U p = 0, (105)
where f d : R + -→ R + is the spectral damping function which is supposed to be increasing and have a sub linear growth at infinity:

f d (λ) ≤ A λ + B, for some A, B > 0. (106)
The operator f d A p is defined in the usual sense of the functional calculus for self-adjoint operator, namely:

∀ n ∈ N, f d A p W n = f d (λ n )W n . (107)
Note that (106) ensures that f d A p is well defined in D(A p ) and that the bilinear form

a p d (U p , U p ) = f d A p U p , U p H ≡ n f d (λ n ) U p , W n H U p , W n H (108)
is well defined in V p = H 1 (ω) 3 while the positivity of f d yields the positivity of the associated quadratic form

∀ U p ∈ H 1 (ω) 3 , a p d (U p , U p ) ≥ 0, (109)
which guarantees the dissipative nature of the additional term in (105). Indeed, the energy identity (100) becomes

d dt E p (U p ) + a p d (∂ t U p , ∂ t U p ) = ω f ∂u p ∂t . (110)
In practice, looking at (107), knowing the {f d (λ n ), n ∈ N} is sufficient, which could be done by matching experimental data. However, it is convenient to use an analytic formula. The following formula is proposed to match experimental data

f d (λ) = α λ + β √ λ + γ, (α, β, γ) ≥ 0 (111)
The reader will note that, as soon as β = 0, f d (A p ) cannot be "identified" to a differential operator and is thus "truly" nonlocal.

Because of this property, it is not easy to rewrite (105) in PDE form. To do so, a solution consists in extending the operator f (A p ) to distributions on ω by duality : for any

U p ∈ D (ω) 3 , we define f (A p )U p ∈ D (ω) 3 by ∀ U p ∈ D(ω) 3 , f (A p )U p , U p = U p , f (A p ) U p (112)
which makes sense because D(ω) 3 ⊂ D(A p ) and is justified by the selfadjoint nature of A p . Modulo this extension, the plate model with dissipation can be rewritten as:

   M p (∇ x ) ∂ 2 t U p + f (A p ) ∂ t U p + A p (x, ∇ x )U p = f ν p , x ∈ ω, t > 0, B p (n, ∇ x )U p = 0, x ∈ ∂ω, t > 0. (113)
Remark 5. For numerical purpose, the dissipation will easily be introduced in the weak formulation (or variational formulation) in space of (105), via the bilinear form a p d . Taking into account (108), this will lead us to use a spectral (or modal) approach for the space discretization of (105) : we shall work in a space spanned by {W n , n ≤ N p } where N p is devoted to tend to +∞.

Inria Modeling a piano.

29

A model for the string-sounboard coupling at the bridge

We shall propose in this section a quite simplified model for the string-soundboard interaction via the bridge. We think that this model is a good starting point to explain the main mechanism of the transmission of string's vibrations to the soundboard.

We shall make the assumptions that the motion of each string, and consequently the one of the attach point at the bridge, remains in the vertical plane containing the string's motion. As a consequence, it is sufficient to give a description of the coupling in a 2D setting. In figure 12, we give a schematic view of the important geometric objects in this plane, denoted (P).

• The segment D p is the intersection of the soundboard ω with the plane (P) : D p = ω ∩ (P),

• The segment line D represents the string at rest. Note that D p and D are not parallel : there is a small angle α between them. This particularity (which is more visible on a violin, for instance, but is also present in a piano) allows the transmission of the longitudinal vibrations of the string to the soundboard.

• The extremity corresponding to the abscissa x = 0 (agraffe) is fixed. The one corresponding to x = L (bridge) is mobile.

• At the bridge extremity (x = L), we suppose a condition of zero torque, which writes (76) in the generalized notations.

• The bridge, assimilated to the (vertical) segment Σ, joins the point x a ∈ ω of the soundboard, that we shall abusively call the "attach point", to the mobile extremity of the string (the real attach point).

• The bridge is supposed to be perfectly rigid (undeformable) and its movement purely vertical along the direction of Σ (as pointed out with the blue arrow on figure 12). In the following, we shall work with the space coordinates (x 1 , x 2 , z) of R 3 chosen in such a way that ω is parallel to the "horizontal" plane (x 1 , x 2 ), the line D p is parallel to x 2 axis and D is parallel to the vertical plane (x 2 , z).

Σ D p D α !"#$%&'()#%*%+),-$%-. v i u i x a
The above assumptions lead to the following kinematic boundary conditions which express that the movement of the string's mobile extremity has no horizontal component while the vertical one is the same than the one of the attach point x a :

   u i (L, t) sin α + v i (L, t) cos α = 0 u i (L, t) cos α -v i (L, t) sin α = u p (x a , t) (114)
The efforts between the two systems string and bridge-soundboard are explained on figure 13, in which the red curve represents the deformed string at a given instant t. More precisely

Chabassier & Chaigne & Joly

• The i th string acts on the soundboard through the string's tension at its extremity T i (L, t),

• The horizontal component of this tension is annihilated by a force of reaction F i r (t) due to the bridge that ensures that the motion of the bridge remains vertical,

• The vertical component of this tension, due to the rigidity of the bridge, is entirely transmitted to the soundboard as a vertical force F i p (t) applied at point x a .

F i r (t) The consequence is the apparition of a right hand side f s = f i s (a point source), for the equation in u p in the soundboard model (see section 4.1 and more precisely (92) and remark 4), :

F i p (t) v i u i x a T i (L, t)
f i s (x, t) = i F i p (t) δ(x -x a ), F i p (t) = cos α T i τ (L, t) + sin α T i (L, t) (115)
where we recall that T i τ and T i are the components of the string's tension in the string's referential, given by (51) . To take into account the fact the bridge is not a 1D object, we shall assume that the force f s is not applied only at a point but distributed on a small portion of the soundboard ω. For this we introduce a positive averaging function χ(x) : R 2 → R + which has a small support (a neighborhood of the origin) and satisfies χ(x) dx = 1. This means that (115) is replaced by

f i s (x, t) = F i p (t) χ(x -x a ), F i p (t) = cos α T i τ (L, t) + sin α T (L, t) (116)
As a consequence, in the purpose of energy conservation, the kinematic conditions (114) have to be modified accordingly. More precisely, the vertical displacement of the mobile extremity of the string will be identified to a mean value of the vertical displacement of the plate ω, using the weight function χ:

       u i (L, t) sin α + v i (L, t) cos α = 0 u i (L, t) cos α -v i (L, t) sin α = ω u p (x, t) χ(x -x a ) dx (117)
In order to rewrite (114) and (116) in a more condensed and abstract form, using the general notation of sections 2.4.2 and 4.1.2, it is useful to introduce the vectors in the space of q-unknowns

τ s α = (sin α, cos α, 0) t ∈ R N , ν s α = (cos α, -sin α, 0) t ∈ R N , (note that ν s α = ν s when α = 0) (118)
in such a way that

F i p (t) ≡ (∇ q H i (∂ x q i , q i ) + Γ ∂ x q i ) • τ s α u i (L, t) sin α + v i (L, t) cos α ≡ q i (L, t) • τ s α , u i (L, t) cos α -v i (L, t) sin α ≡ q i (L, t) • ν s α ( 119
)
Inria Then, the equations for the coupled strings-soundboard model can be rewritten as (we remind that ν p = (1, 0, 0) t ∈ R M , see section 4.1.2 ):

                                         Find q i (t) : [0, L] × R + -→ R N and U p : ω × R + -→ R M , (N = 3, M = 3), such that 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) ν s M p ∂ 2 t U p + f (A p ) ∂ t U p + A p (∇ x )U p = χ(x -x a ) i (∇ q H i (∂ x q i , q i ) + Γ ∂ x q i ) • τ s α ν p 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 + ∇ p H(∂ x q, q) (x b , t) = 0, x b = 0 or L B p (n, ∇ x )U p = 0, on ∂ω, t > 0. (120)
The energy associated to the above strings-sounboard model is naturally defined as

E s,p (q i , U p ) = i E s (q i ) + E p (U p ) (121)
where E s (q i ) is the energy of the i th string (see (59)) and E p (U p ) is the energy of the soundboard (see (99)). The reader will easily verify that any smooth enough solution (q i , U p ) of (120) satisfies the energy identity (see (108) for the definition of a p d (•, •) ):

d dt E s,p (q i , U p ) = - i L 0 R qi • qi - i L 0 Γ ∂ x qi • ∂ x qi -a p d (∂ t U p , ∂ t U p ) (122)
Remark 6. We are conscious that our bridge model is probably oversimplified compared to reality and would deserve to be improved in a future work. For instance, instead of treating the bridge as a part of the soundboard, it would de more realistic to treat it as a (quite rigid) beam sticked to the soundboard and to authorize movements of the bridge which are more complex that the simple vertical translation that we shall authorize here. This would transmit the horizontal transverse vibrations of the string to the soundboard, if the enriched string model of section 2.4.3 were used.

Piano model

By "putting together" the hammer-strings model (87) of section 3 and the strings-sounbard model (120) of section 4.2, one could immediately write a complete model for a piano in the vacuum. The only missing phenomenon to transform this system into a model for a musical instrument is the sound radiation in the surrounding air. This is described by the variations of a pressure field in the air, subject to the standard 3D acoustic wave equation. The pressure field is coupled to the soundboard vibrations through fluid-structure interaction : this is what one calls the structural acoustic equations.

Structural acoustic and sound radiation

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

• ω is, as in section 4.1, the 2D domain in the plane z = 0 that represents the soundboard,

• Ω 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

Ω = R 3 \ Ω f ∩ ω (123)

Chabassier & Chaigne & Joly

We shall distinguish the two sides of the soundboard by introducing ω + ≡ ∂Ω + ∩ ω, ω -≡ ∂Ω -∩ ω, , where Ω ± = Ω ∩ {± z > 0} (124)

In particular, a function q defined in Ω (in H 1 (Ω)) will have two distinct traces on ω, namely q| ω + and q| ω -.

We choose to write the acoustic equations as a first order system and thus introduce as the acoustic unknowns

• the pressure field p(x, z, t) : Ω × R + -→ R,

• the acoustic velocity field V (x, z, t) : Ω × R + -→ R 3 ,

The governing equations are equations describing the interaction of the pressure field (the fluid) with the structure (the soundboard and the furniture part of the piano). The first part of the fluid-structure coupling is given by the kinematic conditions which can be interpreted as boundary conditions for (125) on

∂Ω = ∂Ω f ∪ ω + ∪ ω -
We shall assume that the rim of the piano, Ω f , does not vibrate, which means that we use rigid boundary conditions on ∂Ω f (n f denoting the unit normal vector to ∂Ω f ):

V • n f = 0 on ∂Ω f . (126)
On ω, we express the continuity of the normal velocity of the fluid and the transverse velocity of the plate:

V • e z | ω + = V • e z | ω + = ∂ t u p (≡ ∂ t U p • ν p ) on ∂ω (127)
The second part of the coupling is provided by the description of the efforts applied by the air on the soundboard due to the jump of pressure between the two sides ω + and ω -which constitutes a surfacic force density f a appearing at the right hand side of the plate equation (see (92) and remark 4), namely

f a = p ω = p| ω + -p| ω + . ( 128
)

Inria

Chabassier & Chaigne & Joly counterpart for the coupling with the soundboard through the right hand side. The function χ and the attach point x a are defined in section 4.2. The vectors τ s α and ν s α ∈ R N are defined by (118)) and ν p ∈ R M has been introduced in section 4.1.2.

• Equation (129d), that holds for x ∈ ω and t > 0, governs the movements of the soundboard. For the meaning of the plate unknown U p is given in section 4.1.2. The matrix M p and the operator A p (∇ x ) are given by ( 95) and (96) respectively while the operator f (A p ) is defined in section 4.1.3. The coupling soundboard-strings appears through the first term in the right side, the coupling soundboard-air through the second one.

• Equation (129e), that holds for x ∈ ω and t > 0, corresponds to the boundary condition for the soundboard. The operator B p (∇ x ) is given by (97) and n is the unit normal vector to ∂ω.

• Equations (129f), that holds for (x, z) ∈ Ω and t > 0, governs the variations of the pressure and acoustic velocity fields in the air (see section 5.1).

• Equations (129g) are the interaction equations between the air and the structures, namely with the soundboard and with the rigid part of the piano.

Of course, equations (129) have to be 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) = ∂q i ∂t (x, 0) = 0, x ∈ [0, L],
U p (x, 0) = ∂U p ∂t (x, 0) = 0, x ∈ ω, p(x, z, 0) = 0, V (x, z, 0) = 0, (x, z) ∈ Ω,

We have not studied the mathematical analysis (existence, uniqueness of the solution) of (129, 130) which is probably a hard task and was not our main purpose in this work. However, this system possesses a fundamental stability property through energy estimates. In fact, it is not difficult to establish that any smooth enough solution of (129, 130) satisfies the energy identity (see (108) for the definition of a p d (•, •) ):

d dt E tot (ξ, q i , U p , p, V ) = - i L 0 R qi • qi - i L 0 Γ ∂ x qi • ∂ x qi -a p d (∂ t U p , ∂ t U p ) - i r H Φ H e i ∂ t e i 2 (131)
In (131), the total energy of the system is given by (we remind that Ψ H is a primitive of Φ H , see ( 90))

E tot (ξ, q i , U p , p, V ) = m H 2 ξ 2 + i k H Ψ H e i + i E s (q i ) + E p (U p ) + E a (p, V ). ( 132
)
where E s (q i ) is the energy of the i th string (see (59)), E p (U p ) is the energy of the soundboard (see (99)) while the acoustic energy E a (p, V ) is given by:

E a (p, V ) = Ω ρ a |V | 2 + µ a |p| 2 ) (133)

Conclusion and prospects

The piano model (129) gathers most of the current knowledge concerning the main acoustic and vibratory phenomena occurring in the instrument, but also their reciprocal coupling. When designing this time domain model, a particular effort has been made to keep a strong connection with the fundamental laws of physics, but also geometric and physical parameters. This model couples together the hammer, the nonlinear strings, the soundboard with ribs and bridges, and the radiation of acoustic waves in the free field. However, it is the result of compromises or schematizations of a reality far more complex and subject to specific imperatives that we had sometimes to simplify. Especially,

Inria

Modeling a piano.

35

• Taking the non planar motion of the strings into account would only be relevant if we were able to write a hammer model which explains how his motion is generated, and a bridge model which explains how both polarizations are transmitted to the soundboard. For now, we are lacking experimental data to perform such models, but it seems an interesting extension of our work which should improve the realism of the model.

• The ribs and bridges are considered as heterogeneities of the soundboard's plate model, in terms of material, wood orthotropy direction, but also thickness. However, it is a crude approximation of the reality since the theoretical assumptions of the model are not met in this case (especially, ribs and bridges are not symmetrically arranged on both sides of the soundboard, and the thickness is not varying smoothly). We believe that writing a new plate model that would overcome these assumptions would be a great improvement.

• The hammer is given an initial velocity, which means that we totally neglect the hammer's mechanism and more specifically the flexibility of the hammer's shank. Authors have suggested [START_REF] Askenfelt | From touch to string vibrations. i: Timing in the grand piano action[END_REF] that the pianistic touch lied precisely in this feature, which could seem insignificant at first sight. Therefore, we believe that discriminating different kinds of attacks could be possible if this flexible shank was part of our model.

• Experimental results in [START_REF] Askenfelt | Observations on the transient components of the piano tone[END_REF] show the presence of a sound precursor even before the longitudinal precursor that we mentioned in this article. The authors incriminate the very fast transmission of the key shock waves through the keybed and even the whole piano's cabinet, which radiate in the air, creating the so called "touch precursor". This seems to contribute to the percussive attack of piano sounds.

• It may be trivial to point out that most pianist play several notes at the same time. For now, the model we presented only accounts for the excitation of one note, and assumes that the piano is at rest at the very beginning.

Considering several notes at the same time would be an easy improvement of the model, but would lead to much more heavy computations, which is why we have not presented this feature yet. In the same idea of achieving a better realism, it would then be interesting to model the actions of the pedals and of the dampers.

( a )Figure 1 :

 a1 Figure 1: Schematic view of a grand piano's mechanism.

Figure 2 :

 2 Figure 2: Experimental highlight of frequency dependent and presence of nonlinear phenomena in the piano. Note G3.

  Figure 3(c) represents the spectrum of the longitudinal component of the displacement of the string for such an experiment : the coincidence of peak frequencies between figures 3(b) and 3(c) is enhanced by the presence of the dotted vertical lines, set to integral multiples of 475.0 Hz.

FIG. 15 .FIG

 15 FIG. 12. Time-domain plots of the initial string displacement (upper frame) and sound pressure (lower frame) after striking of the Ao (27.5-Hz) key on the piano. The time scales in the upper and lower frames are synchronized, and the vertical scale units are arbitrary.

Figure 3 :

 3 Figure 3: Experimental highlight of the longitudinal vibrations of the string's contribution to the sound pressure. From [PL88].

Figure 4 :

 4 Figure 4: Spectra of different measured signals when striking the D 1 key, in the range [1 2] kHz.

Figure 5 :

 5 Figure 5: Schematic view of the strings' unknowns.

  (a) Flexural and shear eigenfrequencies of string D 1 Comparison between d'Alembert eigenfrequencies, Timoshenko and its Taylor expansion of formula (19).

Figure 6 :

 6 Figure 6: Eigenfrequencies of a stiff string

  37)from which we infer that the eigenvalues of DF (U ) are the square roots of the eigenvalues of D 2 H(U x ), where D 2 H is the Hessian of H, multiplied by (ρA) -1 . It follows that the local (strict) hyperbolicity of the system (35), namely the Inria Second order expansion H DL2 (u, v) Truncated expansion Happ(u, v)

  65) and (A, B, C) are real N × N matrices with (A, C) symmetric and positive.

Figure 8 :

 8 Figure 8: Schematic view of the additional unknowns.

Chabassier & Chaigne &

  Jolymatrices R and Γ are the diagonal matrices

( a )

 a Geometric description of the hammer (right), based on a real hammer's picture (left). (b) Schematic view of the hammer's crushing e(t) on the string.

Figure 9 :

 9 Figure 9: Schematic description of the hammer.

Figure 10 :

 10 Figure 10: Schematic view of the soundboard's unknowns and geometry.

  (a) Bridges (left) and ribs (right) on both sides of a Fazioli soundboard. (b) Four areas characterizing the presence of bridges and ribs for the Steinway D.

Figure 11 :

 11 Figure 11: Bridges and ribs on soundboards.

Figure 12 :

 12 Figure 12: Schematic description of the bridge (in the plane (P))

Figure 13 :

 13 Figure 13: The forces between the two systems.

Figure 14 :

 14 Figure 14: Geometric configuration of the piano.

RR n°8097

Chabassier & Chaigne & Joly

www.ircam.fr/?&L=1

In the sequel, we will refer to the geometrical configuration of a grand piano, in which the soundboard is horizontal. Inria

Chabassier & Chaigne & Joly truncated version of it, namely U T (t) = χ T (t) U (t) where χ T (t) ∈ [0, 1] is an appropriate cut-off function) namely : f -→ ln | U T (f )|, U T (f ) = +∞ 0

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

Chabassier & Chaigne & JolyInria

The full piano model in PDE form

We build our "complete" piano model by concatenating the content of sections 2.4.2 (or more precisely 2.5), 3, 4.1, 4.2 and 5.1 and more particularly the sub models (87), (120) and (125,126,127). Doing so, we consider a note alone and restrict ourselves the only the strings that are hit by the hammer vibrate. This leads to the following problem:

Strings equations

Strings' boundary conditions -strings / soundboard coupling equations

Soundboard's equations

Soundboard's boundary conditions

Sound propagation

Acoustic boundary conditions -soundboard / air coupling conditions

In the above system, let us recall that N = 3, M = 3 and that

• Equation (129a), that holds for all t > 0, governs the movement and the deformation of the hammer. For the meaning of the unknown ξ(t), the various coefficients, we refer the reader to section 3. The function Φ H is given by (86). The hammer-strings coupling appears through e i and ξ i and x i represents the hammer's impact point on the string.

• Equation (129b), that holds for x ∈ [0, L] and t > 0, governs the movement of the i th string. The meaning of the unknown q i is explained in section 2.4.2 (see (62)). The function H i (p, q) is given, up to the additional index i, by (63) (see also (64)). The matrices M, R and Γ are given by ( 62) and (75) and the function δ H is described in section 3. The hammer-strings coupling appears through F i .

• Equations (129c) are the boundary conditions for the string (we refer to (56) for the definition of J). In particular, the first line represents the kinematic conditions at the point x = L. The second condition is the kinematic