Julien Bensa

Stefan Bilbao

Richard Kronland-Martinet

Julius O Smith

The simulation of piano string vibration: From physical models to finite difference schemes and digital waveguides

A model of transverse piano string vibration, second order in time, which models frequency-dependent loss and dispersion effects is presented here. This model has many desirable properties, in particular that it can be written as a well-posed initial-boundary value problem ͑permitting stable finite difference schemes͒ and that it may be directly related to a digital waveguide model, a digital filter-based algorithm which can be used for musical sound synthesis.

Techniques for the extraction of model parameters from experimental data over the full range of the grand piano are discussed, as is the link between the model parameters and the filter responses in a digital waveguide. Simulations are performed. Finally, the waveguide model is extended to the case of several coupled strings.

I. INTRODUCTION

Several models of transverse wave propagation on a stiff string, of varying degrees of complexity, have appeared in the literature. [START_REF] Hiller | Synthesizing musical sounds by solving the wave equation for vibrating objects[END_REF][START_REF] Bacon | A discrete model of a struck string[END_REF][START_REF] Chaigne | Numerical simulations of struck strings. I. A physical model for a struck string using finite difference methods[END_REF][4][START_REF] Chaigne | On the use of finite differences for musical synthesis. Application to plucked stringed instruments[END_REF] These models, intended for the synthesis of musical tones, are always framed in terms of a partial differential equation ͑PDE͒, or system of PDEs; usually, the simplified starting point for such a model is the one-dimensional wave equation, [START_REF] Fletcher | The Physics of Musical Instruments ͑Springer[END_REF] and the more realistic features, such as dispersion, frequency-dependent loss and nonlinear hammer excitation ͑in the case of the piano string͒, are incorporated through several extra terms. Chaigne and Askenfelt [START_REF] Chaigne | Numerical simulations of struck strings. I. A physical model for a struck string using finite difference methods[END_REF] have proposed such a model ͑see the Appendix for a concise description of this model͒, and used it as the basis for a synthesis technique, through the use of finite differences-the time waveform on a struck piano string is simulated in this way to a remarkable degree of fidelity. 4 Frequencydependent loss, the feature of primary interest in this paper, is modeled through the use of a third-order time derivative perturbation to the dispersive wave equation; while a physical justification for the use of such a term is tenuous, it does give rise to perceptually important variations in damping rates.

In Sec. II, we introduce a model of string vibration, which is of second order in time differentiation; frequencydependent loss is introduced via mixed time-space derivative terms. As it turns out, the model discussed here is a substantial improvement in several different ways. First, the frequency domain analysis of a second-order system is quite straightforward, and it is quite easy to obtain explicit formulas for dispersion and loss curves; this is considerably more complicated for systems which are higher order in time, essentially requiring the factorization of a higher order polynomial dispersion relation. Second, it is easy to prove that our model, when complemented by initial and boundary conditions, is well posed. [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations ͑Wadsworth and Brooks/Cole Advanced Books and Software[END_REF][START_REF] Gustaffson | Time Dependent Problems and Difference Methods ͑Wiley[END_REF] Though we do not give a complete description of this condition here, to say that such an initialboundary value problem is well posed is to say, generally speaking, that solutions may not grow faster than exponen-tially; for linear and shift-invariant systems such as simple stiff string models, the condition can be conveniently expressed in the frequency domain. We show in the Appendix that the PDE first proposed by Ruiz, [START_REF] Ruiz | A technique for simulating the vibrations of strings with a digital computer[END_REF] and later popularized by Chaigne and Askenfelt [START_REF] Chaigne | Numerical simulations of struck strings. I. A physical model for a struck string using finite difference methods[END_REF]4 is in fact not well-posed, and possesses a spurious unstable solution. Third, it becomes easy to develop finite difference schemes, for which precise numerical stability conditions may be easily derived. Finite difference schemes are discussed briefly in Sec. III. Fourth, it is possible to extend the model described here to a more realistic representation of dispersion and loss as a function of frequency through additional terms in the PDE, without compromising well-posedness, or the stability of resulting numerical schemes.

Finally, it is possible to identify the model PDE with a digital waveguide-this filterlike structure models onedimensional wave propagation as purely lossless throughout the length of the string, with loss and dispersion lumped in terminating filters. It thus performs a simulation of a modified physical system.

We show, in Sec. IV A, how one can relate the PDE model presented in Sec. II to a digital waveguide structure, paying particular attention to the relationship between the lumped filters used to model loss and dispersion and the model parameters which define our PDE. In Sec. V, we perform several numerical simulations in order to compare the finite difference and waveguide approaches for this particular problem. In particular, in Sec. VI, we examine in detail a procedure allowing the resynthesis of natural string vibration. Using experimental data obtained from a grand piano, both the terminating filters of a digital waveguide and the parameters of the physical model are estimated over most of the keyboard range. This leads to a simple description of the variation of some of these parameters ͑namely loss parameters and string stiffness͒ over the piano's range which can be used for the convenient synthesis string vibrations at a given excitation point. Finally, in Sec. VII, interstring coupling is discussed and modeled using coupled digital waveguides. This is a further step towards the design of a realistic piano simulator, which should ultimately also model effects such as nonlinear hammer-string coupling and soundboard radiation phenomena.

II. SECOND-ORDER MODELS OF ONE-DIMENSIONAL WAVE PROPAGATION

A. A family of PDEs

Consider a general linear second-order ͑in time͒ wave equation, of the form

ץ 2 y ץt 2 ϩ2 ͚ kϭ0 M q k ץ 2kϩ1 y ץx 2k ץt ϩ ͚ kϭ0 N r k ץ 2k y ץx 2k ϭ0. ͑1͒
Here, y(x,t), the solution, is a function of position xR and time tу0, and q k and r k are real constants; the solution is not uniquely defined until two initial conditions are supplied. Eq. ͑1͒ can be rewritten as

ץ 2 y ץt 2 ϩ2 ͚ kϭ0 M ͑ j␤ ͒ 2k q k ץy ץ t ϩ ͚ kϭ0 N ͑ j␤ ͒ 2k r k y ˆϭ0. ͑2͒
This second-order linear ordinary differential equation with real coefficients will have solutions of the form y ˆ͑␤,t ͒ϭy ˆ0͑␤ ͒e st ͑3͒ for complex frequencies s which satisfy the characteristic polynomial equation

s 2 ϩ2q͑␤ ͒sϩr͑ ␤͒ϭ0 ͑4͒ with q͑␤ ͒ϭ ͚ kϭ0 M ͑ j␤ ͒ 2k q k , r͑␤ ͒ϭ ͚ kϭ0 N ͑ j␤ ͒ 2k r k .
Notice that because only even derivatives appear in the family Eq. ͑1͒, the functions q and r are real. The characteristic polynomial equation has roots

s Ϯ ϭϪqϮͱq 2 Ϫr.
The condition that the initial value problem corresponding to Eq. ͑1͒ be well-posed is that these roots have real parts which are bounded from above as a function of ␤; this is in effect saying that solution growth can be no faster than exponential. A more restrictive ͑and physically relevant͒ condition is that these roots have nonpositive real part for all ␤, so that all exponential solutions are nonincreasing. It is simple to show that this will be true for q͑␤ ͒,r͑ ␤͒у0. ͑5͒

For q and r satisfying Eq. ͑5͒, the imaginary parts of these roots correspond to oscillation frequencies, and the real parts to loss. Clearly, for real wave numbers ␤ such that q 2 рr, the real parts of s Ϯ are simply Ϫq. This case corresponds to normal damped wave propagation; notice in particular that if q depends on ␤, then damping rates will be wave number ͑and thus frequency͒ dependent. If q 2 Ͼr, then both roots are purely real and nonpositive, yielding damped nontraveling solutions. Consider a member of the family defined by Eq. ͑1͒,

ץ 2 y ץt 2 ϭc 2 ץ 2 y ץx 2 Ϫ 2 ץ 4 y ץx 4 Ϫ2b 1 ץy ץt ϩ2b 2 ץ 3 y ץx 2 ץt . ͑6͒
The first term on the right-hand side of the equation, in the absence of the others, gives rise to wavelike motion, with speed c. The second ''ideal bar'' term 10 introduces dispersion, or frequency-dependent wave velocity, and is parametrized by a stiffness coefficient . The third and fourth terms allow for loss, and if b 2 0, decay rates will be frequency dependent. ͓A complete model, for a piano string, is obtained by including a hammer excitation term, f (x,t), possibly accounting for nonlinear effects, on the right-hand side, and by restricting the spatial domain to a finite interval and supplying a realistic set of boundary conditions.͔ This model differs from that of Ruiz 9 in only the last term ͑see the Appendix͒.

The characteristic equation has the form of Eq. ͑4͒, with

q͑␤ ͒ϭb 1 ϩb 2 ␤ 2 , r͑␤ ͒ϭc 2 ␤ 2 ϩ 2 ␤ 4 .
For b 1 , b 2 у0, condition Eq. ͑5͒ is satisfied and this PDE obviously possesses exponentially decaying solutions, and what is more, loss increases as a function of wave number. The PDE of Eq. ͑6͒ possesses traveling wave solutions when q 2 рr, which, for realistic values of the defining parameters for a piano string, includes the overwhelming part of the audio spectrum. For instance, for a C2 piano string, described by the parameters given in Table I, the lower cutoff wave number for traveling waves will be ␤ϭ0.0025, corresponding to a frequency of 0.080 Hz. There is no upper cutoff.

In order to relate this PDE model with a digital waveguide numerical simulation method, it is useful to write the expressions for dispersion and loss directly. Taking s Ϯ ϭϮ j ͑7͒ over the range of ␤ for which traveling wave solutions exist, we obtain

͑␤ ͒ϭϪb 1 Ϫb 2 ␤ 2 , ͑8͒ ͑␤ ͒ϭͱϪ͑ b 1 ϩb 2 ␤ 2 ͒ 2 ϩc 2 ␤ 2 ϩ 2 ␤ 4 .

͑9͒

We will discuss digital waveguide models in detail in Sec. IV A. Finally, we mention that the general form Eq. ͑1͒ serves as a useful point of departure for more accurate models of loss and dispersion in a stiff string. The more terms are included, the better these phenomena can be modeled over the entire frequency range of interest. Though it is difficult to associate physical processes directly with the various extra perturbation terms in the equation, it is at least simple to ensure, through condition Eq. ͑5͒, that the model is wellposed, an important first step in developing stable numerical methods.

B. Boundary conditions

In this section, we provide a brief analysis of pinned boundary conditions, to show that when coupled to our string model, the same analysis of well-posedness may be applied ͑i.e., wave numbers of modal solutions are real͒. Let us now restrict the spatial domain for the problem defined by Eq. ͑6͒ to x͓0,L͔. As Eq. ͑6͒ is of fourth order in the spatial derivatives, we need to supply two boundary conditions at either end, i.e., at xϭ0 and xϭL. Following Chaigne, [START_REF] Chaigne | Numerical simulations of struck strings. I. A physical model for a struck string using finite difference methods[END_REF] we apply pinned boundary conditions,

y͉ xϭ0 ϭy͉ xϭL ϭ ץ 2 y ץx 2 ͯ xϭ0 ϭ ץ 2 y ץx 2 ͯ xϭL ϭ0. ͑10͒
For a solution of the form y(x,t)ϭe stϩ j␤x , from dispersion relation Eq. ͑4͒, there are thus four solutions for ␤ in terms of s,

␤ ϩ 2 ͑ s ͒ϭ Ϫ␥ϩͱ␥ 2 Ϫ4 2 ͑ s 2 ϩ2b 1 s ͒ 2 2 , ͑11a͒ ␤ Ϫ 2 ͑ s ͒ϭ Ϫ␥Ϫͱ␥ 2 Ϫ4 2 ͑ s 2 ϩ2b 1 s ͒ 2 2 ͑11b͒
with ␥ϭc 2 ϩ2b 2 s.

At frequency s, thus, any linear combination y͑x,t ͒ϭe st ͑ a ϩ,ϩ e j␤ ϩ x ϩa ϩ,Ϫ e Ϫ j␤ ϩ x ϩa Ϫ,ϩ e j␤ Ϫ x ϩa Ϫ,Ϫ e Ϫ j␤ Ϫ x ͒ ͑12͒ is a solution to Eq. ͑6͒. Applying the boundary conditions Eq. ͑10͒ to this solution gives the matrix equation

͑13͒

Nontrivial solutions can occur only when det(A)ϭ0, giving the relation

͑ ␤ ϩ 2 Ϫ␤ Ϫ 2 ͒ 2 sin͑␤ ϩ L ͒sin͑ ␤ Ϫ L ͒ϭ0. ͑14͒
Discounting the case ␤ ϩ 2 ϭ␤ Ϫ 2 ͓which yields an identically zero solution y(x,t)], then obviously, solutions are of the form ␤ ϩ ϭn/L, for integer n 0 ͑similarly for ␤ Ϫ ), and the modal frequencies s n are, from the solutions Eq. ͑8͒ and Eq. ͑9͒ of the dispersion relation ͑4͒,

s n ϭ͑n/L ͒ϩ j͑n/L ͒ ͑15͒
over wave numbers for which a traveling solution exists ͑for small b 1 and b 2 , this will be true for all n͒.

III. A FINITE DIFFERENCE SCHEME

In order to solve Eq. ͑6͒ numerically, we may approximate its solution over a grid with spacing X, and with time step T. Equation ͑6͒ can be written as

␦ t 2 yϭc 2 ␦ x 2 yϪ 2 ␦ x 2 ␦ x 2 yϪ2b 1 ␦ t,0 yϩ2b 2 ␦ x 2 ␦ t,Ϫ y ϩO͑T,X 2 ͒, ͑16͒
where the various difference operators are defined by All these operators are ''centered'' about the point (x,t), except for the backward difference operator ␦ t,Ϫ , which is used in order to obtain an explicit algorithm. The approximation is first-order accurate in the time step T, and secondorder accurate in the space step X. ͓It is worth mentioning, that this is but one among many ways of discretizing ͑6͒.͔ We may now rewrite Eq. ͑16͒ as a difference scheme, operating on the grid function y m n , indexed by integer m and n, which will serve as an approximation to y(x,t) at the location xϭmX, tϭnT: Here, the difference scheme coefficients are defined by

␦ x 2 y͑x,
a 10 ϭ͑2Ϫ2 2 Ϫ6 2 Ϫ4b 2 /͒/͑1ϩb 1 T ͒, a 11 ϭ͑ 2 ϩ4 2 ϩ2b 2 /͒/͑1ϩb 1 T ͒, a 12 ϭϪ 2 /͑1ϩb 1 T ͒, ͑18͒ a 20 ϭ͑Ϫ1ϩ4b 2 /ϩb 1 T ͒/͑ 1ϩb 1 T ͒, a 21 ϭ͑Ϫ2b 2 /͒/͑1ϩb 1 T ͒,
where, for brevity, we have introduced the quantities ϭcT/X, ϭT/X 2 .

In order to examine the stability of scheme Eq. ͑17͒, which is, like its generating PDE, Eq. ͑6͒, linear and shift invariant, we may apply frequency domain techniques-in the finite difference setting, these techniques are referred to as Von Neumann analysis. [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations ͑Wadsworth and Brooks/Cole Advanced Books and Software[END_REF][START_REF] Gustaffson | Time Dependent Problems and Difference Methods ͑Wiley[END_REF] This analysis proceeds in a fashion exactly analogous to the analysis applied to the continuous time and space systems of the preceding sections. Shortcutting the process somewhat, we may consider a solution of the form y m n ϭz n e imX␤ , where zϭe sT ͑we could equivalently employ a z transform and a discrete-time Fourier transform͒. We thus obtain the two-step characteristic, or amplification equation

z 2 ϩa 1 ͑ ␤͒zϩa 2 ͑ ␤͒ϭ0,
where the functions a 1 (␤) and a 2 (␤) are defined in terms of the difference scheme parameters of Eq. ͑18͒ by

a 1 ͑ ␤͒ϭϪa 10 Ϫ2a 11 cos͑␤X ͒Ϫ2a 12 cos͑2␤X ͒, a 2 ͑ ␤͒ϭϪa 20 Ϫ2a 21 cos͑␤X ͒.
The necessary and sufficient stability conditions for Eq. ͑17͒ are that the roots of the amplification equation be of magnitude less than or equal to unity, for all wave numbers ␤; for this real-valued quadratic, these conditions can be written in terms of a 1 (␤) and a 2 (␤) as

͉a 1 ͑ ␤͉͒Ϫ1рa 2 ͑ ␤͒р1.
The right inequality is satisfied for b 1 , b 2 у0, and after some algebra, it can be shown that the left inequality is equivalent to the condition 2 ϩ4 2 ϩ4b 2 /р1, further implying that

TрX 2 Ϫ4b 2 ϩͱ16b 2 2 ϩ4͑c 2 X 2 ϩ4 2 ͒ 2͑c 2 X 2 ϩ4 2 ͒ .
The relative ease with which an exact bound such as the above may be derived is a direct consequence of the use of a two-step scheme and the relative simplicity of the model itself; for schemes involving more steps of lookback ͑which results from the discretization of a model with higher time derivatives, such as Ruiz's system͒, this analysis becomes much more involved, though we do approach it nonetheless in the Appendix. This, in addition to reduced memory requirements, is a further advantage of using a second-order model as a starting point.

Equation ͑16͒ is but one of many possible discretizations of Eq. ͑6͒-for instance, replacing ␦ t,Ϫ by ␦ t,0 yields an implicit algorithm, [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations ͑Wadsworth and Brooks/Cole Advanced Books and Software[END_REF] and other implicit schemes such as the -forms discussed in the work of Chaigne 5 may be of interest in reducing numerical dispersion, [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations ͑Wadsworth and Brooks/Cole Advanced Books and Software[END_REF] and may be of higher formal accuracy ͑which may be tempered by the stability re-quirements͒. To emphasize this point, we have plotted the phase velocities and loss curves for the model system of Eq. ͑6͒ versus those of the difference scheme of Eq. ͑17͒ in Fig. 1. Notice that this simple difference scheme is a reasonable approximation to the model only for small ͑i.e., for low frequencies͒. As we will see later in Sec. V, this deviation from the model PDE will account for differences in simulation results obtained from digital waveguide models, which approximates the phase velocity and loss curves directly.

IV. THE DIGITAL WAVEGUIDE MODEL

The digital waveguide approach provides computational models for musical instruments primarily in the string, wind, and brass families. [START_REF] Smith | Digital Waveguide Modeling of Musical Instruments[END_REF] They have also been developed specifically for piano synthesis. [START_REF] Smith | Digital Waveguide Modeling of Musical Instruments[END_REF][START_REF] Smith | Commuted piano synthesis[END_REF][START_REF] Van Duyne | Developments for the commuted piano[END_REF] This section summarizes the basic ideas of the digital waveguide approach, and relates the parameters of a digital waveguide model to an underlying physical model.

A. Background

As mentioned in Sec. I, to arrive at a PDE modeling the piano string, it is fruitful to start with the ideal wave equation and add perturbation terms to give more realistic frequencydependent dispersion and losses. The perturbed PDE is then numerically integrated via a finite-difference scheme ͑or possibly by another approach, such as finite element methods, etc.͒. The digital waveguide approach interchanges the order of these operations: the ideal wave equation is integrated first using a trivial finite-difference scheme, and the resulting solutions are perturbed using digital filters to add frequencydependent loss and dispersion. In the case of a strongly dissipative and dispersive string, the modulus of these so-called loop filters decreases rapidly with frequency, and phase can become strongly nonlinear. For a frequency-domain implementation, this has no effect on computational complexity, but for a time-domain implementation, a larger filter size may be required in order to match the large variations of the phase response.

It has been known since d'Alembert [START_REF] Le Rond D'alembert | Investigation of the curve formed by a vibrating string, 1747[END_REF] that the ideal onedimensional wave equation is solved exactly by arbitrary ͑sufficiently smooth͒ wave shapes propagating in both directions. The digital waveguide formulation works directly with these traveling wave components. An isolated traveling wave is trivially simulated in practice using a delay line. An ideal vibrating string is then modeled as a pair of delay lines, one for each direction of travel.

For digital implementation, the traveling waves are sampled at intervals of T seconds. By Shannon's sampling FIG. 1. Phase velocity ͑top͒ and loss ͑bottom͒ for the model of Eq. ͑6͒ ͑solid line͒ and for difference scheme Eq. ͑17͒ ͑dashed line͒ as a function of the frequency. The model parameters are those corresponding to the note C2, as given in Table I. theorem, [START_REF] Smith | A flexible sampling-rate conversion method[END_REF] the solution remains exact, in principle, at all frequencies up to half the sampling rate. To avoid aliasing, all initial conditions and ongoing excitations must be bandlimited to less than half the sampling rate f s ϭ1/T.

Figure 2 shows the simulation diagram for a digital waveguide model of a rigidly terminated string. The string is excited by the signal E and observed via the signal S. Sampled traveling velocity waves propagate to the right along the upper rail, and to the left along the lower rail. The rigid terminations cause inverting reflections ͑the two Ϫ1 scale factors͒. The delay lines are denoted D i , iϭ1,2,3, and the F i blocks are digital filters to be described further below.

Consider a wavelike solution propagating from a point M 1 to a point M 2 along a string ͑see Fig. 3, top͒. The distance M 1 M 2 will be arbitrarily called l and the propagation time d at the minimal phase velocity. At the observation point M 2 , the wave will have arrived after having undergone the effects of loss and dispersion. In terms of digital waveguides, the wave will undergo a pure delay ͓in the frequency domain, a multiplicative phase factor exp(Ϫid)], times a multiplicative factor F() representing the loss and the dispersion experienced by the wave during this interval ͑see Fig. 3, bottom͒.

Since loss and dispersion are, for this system, linear time-invariant ͑LTI͒ phenomena, even when frequency dependent, the perturbations needed for added realism in the digital waveguide string model are LTI digital filters. Since LTI filters commute, we may lump all of the filtering associated with propagation in one direction into a single LTI filter. These filters are denoted F i in Fig. 2, iϭ1,2,3.

For purposes of computing the output signal S from the input signal E, Fig. 2 may be further simplified.

The two filters labeled F 1 can be replaced by a single filter F 1 2 ͑by commuting one of them with the intervening two delay lines D 1 and Ϫ1 gain͒. A similar simplification is possible for F 3 .

In the same way, the two delay lines labeled D 1 can be replaced by a single delay line D 1 2 ͑having twice the length of D 1 ), and the two D 3 blocks can be replaced by one D 3 2 block.

In general, any uniform section of a linear vibrating string which is excited and observed only at its endpoints can be accurately modeled ͑in one vibrational plane͒ by a pair of digital delay lines, each in series with a digital filter. A discussion of more generalized approaches involving nonuniform string sections, and the relationship with finite difference schemes, is provided in a recent dissertation. [START_REF] Bilbao | Wave and scattering methods for the numerical integration of partial differential equations[END_REF] Since losses and dispersion are relatively weak in vibrating strings and acoustic bores, a low-order filter can approximate very well the distributed filtering ͑infinite-order in prin-ciple͒ associated with a particular stretch of string or bore. ͑We should repeat, however, that the approach outlined above is equally applicable in the case of strongly dispersive or lossy systems, though in these cases, higher-order filters may be necessary.͒ In practice, the desired loss and dispersion filters are normally derived from measurements such as the decay time of overtones in the freely vibrating string. [START_REF] Smith | Techniques for digital filter design and system identification with applicaltion to the violin[END_REF][START_REF] Laroche | Multichannel excitation/filter modeling of percussive sounds with application to the piano[END_REF][START_REF] Va ¨lima ¨ki | Physical modeling of plucked string instruments with application to real-time sound synthesis[END_REF] In the next section, the filter will be derived from the stiff-string model of Eq. ͑16͒. Interestingly, the filter͑s͒ so designed can be mapped back to an equivalent PDE, including many higher-order terms ͑which may or may not have a physical interpretation͒. Lumping of traveling-wave filtering in this way can yield computational savings by orders of magnitude relative to more typical finite difference schemes, [START_REF] Smith | Efficient synthesis of stringed musical instruments[END_REF][START_REF] Karjalainen | Towards high-quality sound synthesis of the guitar and string instruments[END_REF][START_REF] Hikichi | An approach to sound morphing based on physical modeling[END_REF] and this efficiency can be important in applications such as real-time modeling of musical instruments for purposes of automatic sound synthesis or ''virtual acoustic instrument'' performance.

B. Relating digital waveguide parameters to the physical model

We address here the problem of relating digital waveguide filter parameters to the loss and dispersion curves from the physical model discussed in Sec. II. For that, we consider the continuous frequency representation of the loop filter and show its relation with the physical parameters. The digital waveguide parameters can then be obtained by discretization. We do not address here the problem of the time domain implementation of the digital waveguide.

According to Eqs. ͑3͒ and ͑7͒, the transformations of the wave due to propagation along the string segment can be represented in terms of a digital waveguide filter by a multiplicative phase factor exp(sdϩj␤l). Ideally, the modulus and phase of this expression are related to the filter F by ͉F͑ ͉͒ϭ͉e sdϩ j␤l ͉ϭe d ϭe Ϫ͑b 1 ϩb 2 ␤ 2 ͒d , ͑19͒

arg͑F͑ ͒͒ϭarg͑ e sdϩ j␤l ͒ϭdϩ␤l. ͑20͒

In order to write the expressions of the modulus and the phase of the loop filter in terms of the frequency , it is necessary to express the wave number ␤ in terms of . From Eq. ͑9͒, solving for ␤, one gets

␤ 2 ͑ ͒ϭ Ϫ␣ 1 Ϯͱ␣ 1 2 ϩ4␣ 2 ͑ b 1 2 ϩ 2 ͒ 2␣ 2 ͑21͒ with ␣ 1 ϭc 2 Ϫ2b 1 b 2 , ␣ 2 ϭ 2 Ϫb 2 2 .

͑22͒

Since ␤ must be real ͑see Sec. II B͒, we keep only the solution for which the term inside the root is positive. Then,

␤͑͒ϭϮ ͱ Ϫ␣ 1 ϩ ͱ ␣ 1 2 ϩ4␣ 2 ͑ b 1 2 ϩ 2 ͒ 2␣ 2 . ͑23͒
Given that, for realistic piano string modeling, b 1 Ӎ1, b 2 Ӎ10 Ϫ4 , cӍ200, Ӎ1, and Ӎ400, we make the simplifying assumptions

b 1 b 2 Ӷc 2 , b 2 2 Ӷ 2 ͑24͒
͑the second of which also ensures the realness of ␤͒

b 1 2 Ӷ 2 ͑25͒
which permit the following approximation of ␤:

␤͑͒ӍϮ c ͱ 2 ͱϪ1ϩͱ1ϩ4 2 c 4 2 . ͑26͒
In practice, it is helpful to work with more perceptually significant parameters for sound synthesis purposes. For that purpose, we will now suppose that the string is of length L, with perfect reflections at the extremities. The delay D, which corresponds to the propagation time of the minimum phase-velocity wave along the length L can be expressed as

DϭL/cϭ/ 0 , ͑27͒
where 0 is the fundamental frequency ͑rad/s͒ of the ideal string and c is the minimum phase velocity given by cϭ 0 L/. ͑28͒ Thus, D is the propagation time over distance L for a sinusoidal traveling wave tuned to the first resonant mode of the ideal string. This choice of nominal propagation-time D simplifies the frequency-domain computations to follow; however, since phase velocity increases with frequency, the associated propagation filter F will be noncausal in the time domain. This poses no difficulty for frequency-domain implementation.

Next, we may express ␤ as

␤͑͒ӍϮ Lͱ2B ͱ ͑29͒ with ϭϪ1ϩͱ1ϩ4B 2 / 0 2 ͑30͒
in terms of the inharmonicity coefficient 23 B given by

Bϭ 2 0 2 /c 4 , ͑31͒
where c denotes c( 0 ) for notational simplicity.

We now have to choose the sign of ␤ in the expression for the phase. Since we want the output signal to be delayed with respect to the input signal, the loop filter/pure delay combination has to be causal. This means that the phase of the whole transfer function must be negative, i.e.,

ϪDϩarg͑F͑ ͒͒Ͻ0. ͑32͒

This indicates that we choose the negative solution for ␤ in the phase expression. Finally, using Eq. ͑29͒, we arrive at approximate expressions for the modulus and phase of the filter F as a function of frequency ,

͉F͑ ͉͒Ӎexp ͩ ϪD ͫ b 1 ϩ b 2 2 2BL 2 ͬͪ , ͑33a͒ arg͑F͑ ͒͒ӍDϪͱ 2B . ͑33b͒
These expressions serve as the link the PDE model of Eq. ͑6͒ to the lumped filters of the digital waveguide. For the sake of simplicity, one can choose an ''elementary filter''

␦FϭF 1/D such that ͉␦F͑ ͉͒Ӎexp ͩ Ϫ ͫ b 1 ϩ b 2 2 2BL 2 ͬͪ , ͑34a͒ arg͑␦F͑ ͒͒ӍϪ 0 ͱ 2B . ͑34b͒
The filters F of the digital waveguide, which correspond to propagation over a time duration D, can be then easily expressed in terms of ␦F by Fϭ␦F D . ͑35͒

The stability of the digital waveguide model is ensured if the modulus of the filter ␦F is less than one. This condition is here always respected, since the expression in the exponential of Eq. ͑33a͒ is negative.

V. NUMERICAL SIMULATIONS

We now address the validity of the analytical approach of the preceding section to determine digital waveguide loop filters. The modulus and phase of the filters corresponding to the digital waveguide model have been directly linked to the parameters of the physical model. For a given set of PDE parameters, thus, we can design a complete digital waveguide model simulating the signal of the vibrating string at a given location, for a predetermined excitation location. For purposes of comparison, we have generated output signals using both the finite difference scheme discussed in Sec. III, and the digital waveguide model for the same set of model parameters. The digital waveguide has been computed in the frequency domain, allowing the use of Eq. ͑33a͒. The excitation, a Gaussian function of the form e Ϫ͓(xϪx 0 ) 2 / 2 ͔ , simulates an initial velocity distribution of the string at a distance x 0 ϭL/8 from one end of characteristic width ϭK c L H , where K c is an arbitrary constant and L H is the hammer width ͑we do not enter into too much detail here, as hammer modeling is not dealt with in this paper͒. The signal is observed at the location 9L/10. We have performed simulations for the notes C2, C4, and C7 using the parameters in Table I ͑all of which are taken from values provided in the papers by Chaigne and Askenfelt, [START_REF] Chaigne | Numerical simulations of struck strings. I. A physical model for a struck string using finite difference methods[END_REF]4 except for the parameter b 2 , which comes from the calibration procedure applied subsequently in Sec. VI.

Figure 4 shows the two first periods of the waveforms generated by both approaches. The amplitudes are similar. Nevertheless, there is a slight discrepancy between the two signals, due to the numerical dispersion introduced by the finite difference scheme ͑see Sec. III͒. This discrepancy can be better seen by comparing the phase velocity of the two signals in Fig. 5. This corresponds to the phase velocity plotted in Fig. 1. The phase velocity of the signal produced by the waveguide is similar to the one of the model.

The long-time behavior of the generated signal is also very similar. In Fig. 6 spectrograms obtained over the whole length of the sound are shown. It is clear that the global damping behavior is similar in the low-frequency range. However, at high frequencies, the finite-difference model suffers from an artificially high propagation gain, as derived in Sec. III. The fundamental frequencies are essentially equal, but the wave dispersion due to string stiffness is different due to the numerical dispersion introduced by the difference scheme. ͑We can see in Fig. 7 a slight tuning difference in the high-frequency partials.͒ In summary, these figures illustrate the extent to which the waveguide model provides a more accurate digital simulation of stiff, lossy strings with respect to both attenuation and dispersion of wave propagation, when compared with finite difference schemes.

VI. CALIBRATION OF PHYSICAL MODEL PARAMETERS FOR A GRAND PIANO FROM EXPERIMENTAL DATA

A. Experimental setup

In order to calibrate the parameters of the physical model, data were collected on an experimental setup consisting of a Yamaha Disklavier C6 grand piano equipped with sensors ͑see Fig. 8͒. The string vibrations were measured using an accelerometer located at the bridge level. The Disklavier allows the piano to be played under computer control. These measurements were made at the Laboratoire de Me ´canique et d'Acoustique, Marseille, France.

In order to ensure that the measurements were taken under similar excitation conditions, we measured the hammer velocity using a photonic sensor ͑MTI 2000, probe module 2125H͒ ͑see Fig. 9͒.

Since we were interested in exciting a large portion of the frequency spectrum while remaining in the linear domain of vibration, we chose a medium ͑mezzo-forte͒ hammer velocity of 2.2 m s Ϫ1 , which corresponds roughly to a MIDI value of 80. Such a hammer velocity allows the generation of about 140 spectral components for low frequency tones with a reasonable signal-to-noise ratio ͑see Fig. 10͒.

For each note of the grand piano, the optical sensor was placed close to the hammer and the accelerometer ͑B&K 4374͒ at bridge level. For notes corresponding to double or triple sets of strings, the accelerometer was placed as close as possible to the strings ͑see Fig. 11͒. Due to the imprecision of MIDI coding, several measurements were made, until a target value of the hammer speed was obtained. We have deemed an uncertainty of Ϯ0.1 m s Ϫ1 for the hammer velocity to be acceptable, as the estimation of modal frequency and decay rates is relatively insensitive to such an error. The acceleration measured at the bridge level was then digitally recorded at 16 bits at a sampling rate of 44.1 kHz, before being entered in the database.

B. Estimation of parameters

Because this model is intended for use in the context of musical sound synthesis, we here discuss the calibration of b 1 and b 2 , and the determination of the stiffness parameter for a given string.

FIG. 8. Experimental setup. The grand piano was isolated in an anechoic room, and both the string vibrations and the acoustical radiated signal were measured. The string vibrations were measured using both an accelerometer located at the bridge level and a laser vibrometer, while the acoustic signal was measured at the head level of the pianist using an artificial headset. Our library of measured data also includes signals corresponding to various hammer velocities ͓referred to in recent and forthcoming articles concerned with the hammer-string interaction ͑Ref. 31͔͒, but for this paper, we only need acceleration measurements for each string. To estimate the damping factor associated with each modal component of the signal, we used a signal processing technique based on the theory of analytic signals. [START_REF] Gabor | Theory of communication[END_REF][START_REF] Kronland-Martinet | Modelling of natural sounds using time-frequency and wavelet representations[END_REF] The analytic signal representation provides an easy way of extracting both the instantaneous frequency and the damping of each modal component, through band-pass filtering. To isolate each component in frequency, we used a truncated Gaussian window, the frequency bandwidth of which was chosen so as to minimize smoothing effects over the attack duration and to avoid overlapping two successive frequency components. The Gaussian was employed since its timebandwidth product is minimized. As a consequence, it optimizes the exponential damping support after convolution with a causal single component for a given band pass filtering. [START_REF] Aramaki | Resynthesis of coupled piano string vibrations based on physical modeling[END_REF] The analytic signal associated with each individual component facilitates the estimation of both the instantaneous frequency and the amplitude modulation law of the component. The frequency dependent damping factor is directly related to the amplitude modulation law of each partial. According to the physics of a single string vibrating in one plane, the amplitude modulation of each component is expected to be exponential. This makes the damping factor easy to estimate by the measurement of the slope of the logarithmic representation of the amplitude of the analytic signal. [START_REF] Va ¨lima ¨ki | Physical modeling of plucked string instruments with application to real-time sound synthesis[END_REF] This technique has advantages over other parametric methods such as Prony's method, [START_REF] Markel | Linear Prediction of Speech ͑Springer[END_REF] mainly due to its ability to extract a coherent mean damping factor when multiple components are present. In fact, the hammer usually strikes not one, but two or three strings simultaneously. The coupling gives rise to perceptually significant phenomena such as beating and two-stage decay; 28 these effects are not accounted for in our model Eq. ͑6͒. For these multistring notes, the calculated damping coefficients can be thought of describing the global perceived decay of the sound.

For the same reason, the spectral representation of each partial is the result of the summation of several contributions due to each string and each polarization of vibration. The phase of the analytic signal allows an accurate estimation of the mean frequency of the partials. Actually, it permits the calculation of the instantaneous frequency which is a timedependent function oscillating around a mean value. This mean value coincides with the spectral centroid of the partial [START_REF] Ph | Characterization of acoustic signals through continuous linear time-frequency representations[END_REF] and consequently with the more likely perceived frequency. From the mean frequency values estimated this way for each partial, B may be deduced for each note, and consequently 0 , the fundamental frequency of the corresponding ideal string. The inharmonicity factor B is plotted as a function of the frequency in Fig. 12; B is an increasing function of the note number, except over the bass range, where the strings are double-wrapped ͑this behavior has also been measured by Conklin 30 ͒.

Combining Eqs. ͑8͒ and ͑29͒, one obtains

͑ ͒ϭϪb 1 Ϫb 2 ͩ 2 2BL 2 ͓Ϫ1ϩ ͱ1ϩ4B ͑ / 0 ͒ 2 ͔ ͪ .

͑36͒

Then, from the value of obtained for each partial, b 1 and b 2 may be estimated for a given tone. The evolution of these parameters as a function of the fundamental frequency is shown in Fig. 13. We see that b 1 and b 2 are both increasing functions of MIDI note number, indicating increasing loss as one approaches the treble range. In Fig. 13, we have also fit extremely simple curves to the loss parameter data. The fits are linear as a function of the fundamental frequency, and are given by

b 1 ϭ4.4ϫ10 Ϫ3 f 0 Ϫ4ϫ10 Ϫ2 , ͑37a͒ b 2 ϭ1.0ϫ10 Ϫ6 f 0 ϩ1ϫ10 Ϫ5 . ͑37b͒
These simple empirical descriptions of b 1 and b 2 allow the reproduction of piano tones whose damping will be close to that of the perceived acoustic note. A multistring waveguide model has also been designed, [START_REF] Aramaki | Resynthesis of coupled piano string vibrations based on physical modeling[END_REF][START_REF] Bensa | Note and hammer velocity dependence of a piano string model based on coupled digital waveguides[END_REF] allowing for beating and two-stage decay. Its use for synthesis purposes is discussed in Sec. VII. The deviations in the curves ͑Fig. 13͒ from the linear fits must be attributed to the impedance ratio between the strings and the soundboard, which varies over the length of the bridge. As a result, the curves should not be interpreted as impedance curves as they do not represent measurements at a single point and data are taken only at string modal frequencies. For a detailed discussion of soundboard impedance and its measurement, we refer to the work of Giordano. [START_REF] Giordano | Mechanical impedance of a piano soundboard[END_REF][START_REF] Giordano | Sound production by a vibrating piano soundboard: Experiment[END_REF]

VII. WAVEGUIDES AND SOUND SYNTHESIS

The determination of b 1 , b 2 , B, and 0 for each note allows for an explicit expression of the behavior of the filter F as a function of note number. In order to represent the evolution of the elementary loop filter in terms of note value, we show in Fig. 14 the modulus and phase of the elementary filter ␦F, normalized with respect to the time delay D, ␦FϭF 1/D .

͑38͒

In order to understand the general behavior of the modulus and the phase of the loop filter with respect to the note played, we expand their expressions ͓Eq. ͑34a͔͒ for 4B(/ 0 ) to third order near zero. We obtain

͉␦F͑ ͉͒Ӎexp ͩ Ϫ ͫ b 1 ϩb 2 2 c 2 ͬͪ , ͑39a͒ arg͑␦F ͒Ӎ 3 B 2 0 2 . ͑39b͒
The modulus ͑which also accounts for the losses at the endpoints͒, is decreasing with note number as shown in Fig. 14. This is mainly due to the increase in b 1 . But the param-increasing, leading to a modification of the slope of the modulus versus frequency for different note numbers. We note that this behavior is slightly different for the wrapped strings ͑A0, A1͒ than for the other strings ͑A2, A3͒. We also note that although B is mainly an increasing function of note number, the phase of the filters ␦F grows less rapidly for the bass tones than for higher tones. This is due to the fact that the phase of the filter depends not only on B, but also on the fundamental frequency, as evidenced by Eq. ͑39b͒. Though this expression is meaningful only for the first few partials, it is clear that phase dispersion decreases as a function of fundamental frequency.

In the case of the piano, the string is struck at a distance of approximately one-eighth to one-sixteenth of its length from the bridge, depending on the note. We are only interested in the vibration generated at the bridge termination, as this is the mechanism by which energy is transmitted to the soundboard. This situation corresponds to a digital waveguide structure identical to the one presented in Fig. 2, except that the loop filters F 2 and F 3 are combined. Using the parameters estimated by experiment, one can reproduce the vibration generated by a single string. Figure 15 shows the evolution of the amplitude of the first six partials of the vibration velocity for the note E1, respectively, measured on the piano and generated by the digital waveguide model. From a perceptual point of view, the two sounds are identical.

If the tones are produced by two or three strings struck simultaneously, the basic digital waveguide model still generates a signal having the same amplitude and damping of the modes. It does not, however, account for the modulations and double decays due to the coupling of the strings at the bridge. [START_REF] Weinreich | Coupled piano strings[END_REF] Using two or three coupled digital waveguides, and thus allowing for energy transfer between the strings, one can easily reproduce this phenomenon. [START_REF] Smith | Efficient synthesis of stringed musical instruments[END_REF][START_REF] Aramaki | Resynthesis of coupled piano string vibrations based on physical modeling[END_REF][START_REF] Bensa | Note and hammer velocity dependence of a piano string model based on coupled digital waveguides[END_REF] Figure 16 shows the time evolution of the amplitude of the first six partials for the tone C2 ͑two strings͒, using the coupled- eter b 2 , which allows for frequency dependent loss is also strings model described in a recent publication. [START_REF] Aramaki | Resynthesis of coupled piano string vibrations based on physical modeling[END_REF] The modulations are essentially perfectly reproduced. Moreover, this model follows directly from the physics of coupled strings. In fact, the loop filters are again related to the parameters of the physical model, and numerous sound transformations are conceivable. One could, for instance, extend the use of such a waveguide model to practically unrealizable situations involving, for example, widely mistuned strings or coupled strings with differing material properties.

VIII. CONCLUSIONS

We have presented a model of transverse vibrations on a string which includes the effects of stiffness and frequencydependent loss. This model possesses several advantages over those proposed previously, in particular that it can be framed as a well-posed initial-boundary value problem ͑leading to stable numerical methods͒, and also that it can be easily related to digital waveguides. The source of these good properties is the fact that this model allows only two traveling-wave type solutions; the nonphysical third unstable term in the model of Ruiz ͑which can lead to difficulties both analytically and numerically͒ is thus eliminated in favor of higher-order spatial terms. It is also simple to write expressions for dispersion and loss as a function of frequency in terms of the model system parameters-such information is critical for the design of the terminating filters in a digital waveguide implementation.

For the sake of comparison, we have performed numerical simulations of the model system, with pianolike parameters, using both finite differences and a digital waveguide; the most notable distinction is the complete lack of numerical dispersion ͑which leads to mode mistuning͒ in the waveguide implementation. On the other hand, the waveguide allows the computation of a solution ͑''sound''͒ only at preselected points on the string, whereas a finite difference scheme computes the entire string state in ͑sampled͒ physical form. This is not a drawback for sound synthesis applications because, only the behavior of the string at the bridge is of interest in most stringed instruments. Moreover, physically accurate outputs from additional points along the string are easily added to a digital-waveguide simulation at the price of one small digital filter each.

A set of experimental data measured from a grand piano was used in order to calibrate the PDE model parameters over the entire keyboard range. String vibration was measured at the bridge through the use of an accelerometer, for each note on the piano, and for an average hammer velocity. The piano employed was equipped with sensors to provide hammer velocity data; from these measurements, all the parameters relating to the relevant PDE model were estimated. Given that the model itself is not completely physical-that is, the various loss mechanisms, interstring coupling, as well as energy transfer to the soundboard are modelled, for simplicity, as internal to the string itself-these parameters must be considered as those describing an ''equivalent'' string, under fixed termination. The equivalent parameters, are, however, sufficient for the resynthesis of piano tones to a high degree of fidelity, when a digital waveguide is employed. The digital waveguide model was also extended in order to directly take into account the effects of interstring coupling, through the use of two or three coupled waveguides.

The modeling of the excitation mechanism for the piano string ͑i.e., the hammer͒ is also of great importance, and must be carried out with some care; we have not addressed this issue here. As has been shown in the work of Chaigne and Askenfelt, it is possible to design a nonlinear hammer, which, when applied to a stiff string with frequencydependent loss, produces signals quite similar to those measured on a real piano. The problem of extracting hammer parameters from measured data is also worthy of future research.

little time here explaining the significance of the difference, which has a radical effect on the analysis of the system as a whole.

We can examine the well posedness of the system by inserting a solution of the form y(x,t)ϭe stϩ j␤x into Eq. ͑A1͒, in order to obtain a dispersion relation,

Ϫ2b 3 s 3 ϩs 2 ϩ2b 1 sϩc 2 ␤ 2 ϩ 2 ␤ 4 ϭ0. ͑A2͒
This is a cubic in the variable s ͓the quantity on the left-hand side is referred to as the symbol 7 of Eq. ͑A1͔͒, and again, as discussed in Sec. II, a necessary condition for the system of Eq. ͑A1͒ to be well posed is that the roots of this equation have real part bounded from above as a function of spatial frequency ␤. It is simple to see that the real part of at least one root of Eq. ͑A2͒ will be positive and unbounded as a function of wave number. As this is a third-degree polynomial equation with real coefficients, one root will always be real, and the two others occur as a complex conjugate pair ͑or perhaps as two other real roots͒. Consider Eq. ͑A2͒ as ͉␤͉ becomes large. In this case, the three roots will behave as

Ϫ2b 3 s 3 ϩ 2 ␤ 4 Ϸ0
and will be evenly distributed over a circle of radius ( 2 ␤ 4 /(2b 3 )) 1/3 . If b 3 Ͼ0 ͑as suggested in the numerical experiments in the papers by Chaigne and Askenfelt 4 ͒, then there will be one positive real root of the magnitude mentioned above, clearly unbounded as a function of wave number ␤. ͑If b 3 Ͻ0, there will two roots in the right half plane, of this same magnitude, at an angle approaching Ϯ60 degrees with respect to the positive real axis.͒ We have thus shown that the initial value problem corresponding to the system of Eq. ͑A1͒ is, formally speaking, ill-posed. We can extract some more detailed information by asking under what conditions the roots of Eq. ͑A2͒ have positive real part. A straightforward application of the Routh-Hurwitz stability criterion [START_REF] Van Valkenburg | Network Analysis ͑Dover[END_REF] to Eq. ͑A2͒ shows that in fact, if b 3 Ͼ0, there is always exactly one real positive root, regardless of the value of the wave number ␤.

The following question then arises: How can we explain the apparently stable behavior exhibited by simulations 4 of these equations? Indeed, for realistic piano string parameters, the numerical integration routine provided in the paper by Chaigne and Askenfelt 4 is stable, and produces piano sounds of excellent quality. A first guess might be that the above analysis is incomplete due to the neglect of boundary conditions. Using the boundary conditions supplied by Chaigne, however, leads to an analysis identical to that performed in Sec. II B-the modal frequencies for the string system defined by Eq. ͑A1͒ will be given by solutions of the dispersion relation Eq. ͑A2͒ under the replacement of ␤ by n/L for integer n. For any n, there will be exactly one modal frequency s n with positive real part. Thus the instability persists even in the presence of boundary conditions.

We must then conclude that discretization has a stabilizing effect on system of Eq. ͑A1͒. To explore this idea in more detail, consider the discretization, 4 which can be written as The solution to the recursion will be bounded and decay if the solutions to this equation are confined to the interior of the unit circle for all ␤͓Ϫ/X,/X͔. It is simple to show that this is in fact true, for any of the choices of parameters given in the papers by Chaigne and Askenfelt. 4 This does not mean, however, that the difference scheme can be considered to be numerically stable in the Von Neumann sense. [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations ͑Wadsworth and Brooks/Cole Advanced Books and Software[END_REF] This is a rather subtle point, and is worth elaborating. According to the Lax-Richtmeyer equivalence theorem, [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations ͑Wadsworth and Brooks/Cole Advanced Books and Software[END_REF] if the initial-boundary value problem is wellposed, the solution to a finite difference scheme will converge to the solution of the model problem if it is consistent and stable. In this case, though, the model system is not well posed, and thus no finite difference can possibly converge to a stable solution in some limit as the time step T and the grid spacing X approach zero. The difference scheme Eq. ͑A3͒ is indeed consistent with ͑A1͒ to first-order accuracy ͑and we note that if one does wish to use this ill-posed model system, it is in fact possible to design second-order accurate explicit methods͒, but it is possible to show ͑as we expect͒ that it cannot be stable in the limit as T becomes small. Because the recursion is of third order, the analysis is somewhat involved, and requires the application of the Schur-Cohn recursive procedure [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations ͑Wadsworth and Brooks/Cole Advanced Books and Software[END_REF] ͑the discrete time analog of the Routh-Hurwitz stability test, again allowing us to check the stability of a polynomial without explicitly calculating the roots͒. Nevertheless, it is possible to show in this way that a necessary condition that the roots of the polynomial of Eq. ͑A5͒ be inside the unit circle is that b 3 /Tр 2 ϩ4 2 р1. ͑The second inequality is exactly the necessary stability condition given in the paper by Chaigne and Askenfelt. [START_REF] Chaigne | Numerical simulations of struck strings. I. A physical model for a struck string using finite difference methods[END_REF] ͒ Al-though we have not provided all the details, we note that it suffices to perform the Schur-Cohn test at the wave number ␤ϭ/X in order to arrive at these conditions. Clearly, these two conditions cannot be satisfied if Tрb 3 ⇒Sample rateу1/b 3 and thus for a small enough time step, the system poles must cross to the exterior of the unit circle, regardless of the grid spacing X.

It should be said, however, that because b 3 is in general quite small for realistic piano string models ͑on the order of 10 Ϫ9 ), for any reasonable sample rate in the audio range, the recursion does not exhibit this unbounded growth. On the other hand, as we have shown in this paper, it is simple enough to dispense with the nonphysical solution and all the concomitant analysis by making use of a simpler secondorder model.

y m nϩ1 ϭa 10 y m n ϩa 11

 11

FIG. 2 .FIG. 3 .

 23 FIG. 2. Digital waveguide model of a rigidly terminated string.

FIG. 6 .FIG. 4 .

 64 FIG.6. Spectrograms of the output of the finite difference scheme ͑at top͒ and the waveguide model ͑at bottom͒ for the note C2.

FIG. 5 .

 5 FIG.5. Phase velocity for the waveguide model ͑solid line͒ and for difference scheme as a function of the frequency for the note C2. The analytically obtained phase is not shown, as it is identical ͑by definition͒ to the phase response of the waveguide network.

FIG. 9 .

 9 FIG.9. Optical sensor used to measure the hammer velocity. The velocity is obtained through the duration corresponding to the travel time of the hammer between two reflectors placed on it.

FIG. 12 .

 12 FIG. 12. The measured inharmonicity factor B.

FIG. 14 .

 14 FIG. 14. Normalized modulus and phase of F for selected tones.

y m nϩ1 ϭa 10 y m n ϩa 20 y m nϪ1 ϩa 11

 11 involves three steps of lookback, reflecting the degree of the model system of Eq. ͑A1͒. Here, the difference scheme coefficients are defined bya 10 ϭ͑2Ϫ2 2 Ϫ6 2 ϩb 3 /T ͒/D,a 20 ϭ͑Ϫ1ϩb 1 Tϩ2b 3 /T ͒/D, a 11 ϭ͑ 2 ϩ4 2 ͒/D, ͑A4͒ a 12 ϭ͑b 3 /TϪ 2 ͒/D, a 21 ϭa 30 ϭ͑Ϫb 3 /T ͒/D, where, for brevity, we have again used ϭcT/X, ϭT/X 2 , as well as Dϭ1ϩb 1 Tϩ2b 3 /T.Let us now examine the characteristic polynomial, which can be written asz 3 ϩa 1 ͑ ␤͒z 2 ϩa 2 ͑ ␤͒zϩa 3 ͑ ␤͒ϭ0 ͑A5͒with a 1 ͑ ␤͒ϭϪa 10 Ϫ2a 11 cos͑␤X ͒Ϫ2a 12 cos͑2␤X ͒, a 2 ͑ ␤͒ϭϪa 20 Ϫ2a 21 cos͑␤X ͒, a 3 ͑ ␤͒ϭϪa 30 .

  ͑For the moment, we concentrate on the pure initial value problem and ignore boundary conditions-we will return to this subject in Sec. II B.

͒ Because this equation describes a linear and shift-invariant system, it is possible to analyze it through Fourier techniques. Defining the spatial Fourier transform y ˆ(␤,t) of y(x,t) by y ˆ͑␤,t ͒ϭ 1 ͱ2 ͵ Ϫϱ ϱ y͑x,t ͒e Ϫ j␤x dx,

TABLE I

 I

		. Physical model parameters for piano tones C2, C4, and C7.
		C2	C4	C7	Units
	L	1.23	0.63	0.10	m
	c	160.9	329.6	418.6	m s Ϫ1
		0.58	1.25	1.24	m 2 s Ϫ1
	b 1 b 2 F s	0.25 7.5ϫ10 Ϫ5 16 000	1.1 2.7ϫ10 Ϫ4 32 000	9.17 2.1ϫ10 Ϫ3 m 2 s Ϫ1 s Ϫ1 96 000 s Ϫ1

ACKNOWLEDGMENTS

This work was supported in part by the DSP ͑France͒ under Contract No. 016060 and by a fellowship from Region PACA ͑France͒.

APPENDIX: THE PIANO STRING MODEL OF CHAIGNE AND ASKENFELT

The results in this section have appeared, in a similar form, in the thesis of Ruiz. [START_REF] Ruiz | A technique for simulating the vibrations of strings with a digital computer[END_REF] We have added various comments regarding well-posedness and numerical stability.

The stiff string model in the thesis of Ruiz [START_REF] Ruiz | A technique for simulating the vibrations of strings with a digital computer[END_REF] and in the papers by Chaigne and Askenfelt 3,5 is described by the following equation:

͑A1͒

This model differs from Eq. ͑6͒ only by the replacement of the term 2b 2 ץ( 3 y/ץx 2 ץt) by 2b 3 ץ( 3 y/ץt 3 ); this model also allows for frequency-dependent loss, but the system itself is of a quite different character, due to the increased degree of the equation with respect to the time variable. We spend a FIG. 16. Amplitude of partials one to six for the note C2 as a function of time and frequency.


Didier Marot

Alain Alexis

Pierre Thomas

TOWARDS A BETTER FORECAST FOR RESERVOIR SLOPE SLIDING

Keywords: Soft soils, Mud, Bank, Sliding, Dam reservoir, Emptying

The siltation of reservoir causes many problems of management and environmental disturbances that are important to minimize. Solving these problems requires in particular a good knowledge of mud bank instabilities. In partnership with the Laboratoire National d'Hydraulique (LNH), a subsidiary of Electricité de France (E.D.F.), authors realize a study of soft soils sliding entailed by reservoir emptying. Using a geotechnical flat sliding mechanism, the theoretical formulation of the safety factor can be established when the sediment mass is submerged and emergent. The measurement campaign carried out makes it possible to obtain bathymetric and topographic surveys and a geotechnical characterization of the soil for the whole reservoir of Grangent dam, located on the Loire (France). The proposed model appears to agree satisfactorily with the instability of the zones emerged because of the emptying and with the stability of the profile after emptying. The use of this model within the framework of forecasts makes it possible to calculate significant minimal displaced heights.

INTRODUCTION

Hydroelectric dams and their storage reservoirs allow for the production of electricity, the creation of water reserves for consumption or irrigation, and the regulation of river flow rate. The long-term future of these installations requires the dredging (often non-viable economically) or the emptying of storage reservoirs. This release of stored water is carried out with two objectives:

-periodic inspection of the structure and the carrying out of maintenance work, -expulsion (flushing) of sediments accumulated as upstream slopes of the dam. The improvement of the sedimentary control and the development of operational solutions require increased knowledge of the impacts of emptying a storage reservoir on the slope stability of its silt deposition. [START_REF] Bouchard | Mud investigations and management simulation of reservoir sediments[END_REF], [START_REF] Paquier | River bed deformation calculated from boundary shear stress[END_REF] presented numerical codes which simulate erosion and deposition in a storage reservoir at times of high waters. However, emptying is also accompanied by a considerable amount of slide of the upstream slopes. [START_REF] Zhou | One-dimensional mathematical model for suspensed sediment by lateral integration[END_REF] arbitrarily set a maximum slope gradient for the sediments, but the difficulty lies in the evaluation of this maximum slope as this is influenced by several factors:

-the thickness of the sediment mass -the slope of the bedrock -the soils characteristics The objective of this study is thus to propose a criteria for upstream slope slide during emptying of storage reservoirs of dams, and to then validate it on a real site: the Grangent dam (Loire-France). An initial criteria is developed with a rectangular-shaped slice [START_REF] Alexis | Mud bank sliding in a dam reservoir[END_REF]. It supposes that the shear plane and the water-sediment interface are parallel. In order to deal with a case more typical of the natural configuration while freeing itself from this simplistic assumption, we develop a sliding criteria for a slice of sediment with a trapezoidal shape. Moreover to improve the validation of proposed criteria, the topographic profiles carried out before the construction of the dam are exploited.

BIBLIOGRAPHICAL ASPECTS

Landslides can be classified into two great rheological families [START_REF] Filliat | La pratique des sols et fondations[END_REF]: sliding in the traditional sense of the word and mud flows. The move to the second type may in particular be obtained by increases in the water content. We therefore need to model the falling of upstream slopes in accordance with these two approaches.

Hydraulic Flow Criteria (Mud Flows)

The low consistency of the mud makes it possible to consider it as a fluid. Its behaviour can be predicted by various models, including viscoelastic models [START_REF] Feng | Mudbanks of the southwest coast of India IV: mud viscoelastic properties[END_REF] and viscoplastic models [START_REF] Coussot | On the behaviour of fine mud suspensions[END_REF]. The plasticity of material results in the need to reach a threshold stress value in order to start the movement. [START_REF] Coussot | Spécificités de l'hydraulique en torrent de montagne[END_REF] specified that in order for flow to take place, it is necessary that the thickness of the fluid h e should satisfy the following inequality:

α sin γ τ h sat s e ≥ (1)
τ s : threshold shear stress; γ sat : bulk density of the saturated fluid; α: inclination of the slope.

This model supposes that the shear plane and the water-sediment interface are parallel.

When the shearing of the fluid gets underway, the velocity profile depends on the behavioural law selected [START_REF] Suhayda | Explanation of submarine landslide morphology by stability analysis and rheological models[END_REF][START_REF] Périgaud | Mécanique de l'érosion des silts[END_REF][START_REF] De Kee | Flow of viscoplastic fluids on an inclined plane: evaluation of yield stress[END_REF].

Geotechnical Slip Criteria

Although of low consistency, the mud can also be regarded as a soil. The study of the stability of the sediment massif is then carried out according to two approaches: the finite element method, which is a numerical method of calculation of the variational type [START_REF] Zou | Search for critical slip surfaces based on finite element method[END_REF], and the limit failure state design. The limit failure state design involves studying the stability of the sediment mass at the boundary equilibrium, i.e. when the failure is on the point of occurring. This study is carried out by the introduction of the safety factor (F) which is the ratio of the shear strength to the shear stress. The unstable volume can be divided into a certain number of vertical slices and the boundary equilibrium of each slice is studied. Two types of calculation are then possible. The first is based on the formulation of the safety factor for an arbitrarily defined sliding surface. This study can be undertaken using two dimensional models [START_REF] Okagbue | Engineering properties and embankments stability of road cuttings along the Enugu-Onitsha express road, southeastern Nigeria[END_REF] and can also take account of the evolution of the pore pressure which influences the stability of the sediment mass [START_REF] Fredlund | Analyses for the stability of potash tailings piles[END_REF]. [START_REF] Lam | A general limit equilibrium model for three-dimensional embankments stability analysis[END_REF] have for their part carried out a three-dimensional analysis. The second type of calculation by the slice method is based on researching the critical sliding surface, obtained by minimizing the safety factor. This safety factor and the corresponding critical sliding slope gradient can then be found by various numerical methods, in two dimensions [START_REF] Chen | Random trials used in determining global minimum factors of safety of embankments[END_REF] or even in three dimensions [START_REF] Yamagami | A search for the critical slip surface in three-dimensional embankments stability analysis[END_REF].

The two dimensional general expression of the equilibrium of each slice leads to four equations: . -2 equations of force equilibrium, -1 equation of moment equilibrium, -1 equation of the Coulomb criterion. Thus for a sediment mass which comprises n slices, 4n equations are obtained (the three dimensional study results from a generalisation of the equations obtained). The resolution of the equilibrium of all the slices leads to (6n -2) unknown factors. There are thus (2n -2) equations missing to solve the problem, and we must introduce various simplifying assumptions.

An assumption is generally made: the normal force applied to the lower surface of the slice is assumed to apply to the centre of this surface. The number of unknown factors is thus reduced to (5n -2). Consequently (n -2) assumptions are necessary to the determination of all unknown factors. These assumptions generally relate to the shape of the sliding surface, the application point of the inter-slice forces or to the expression of these inter-slice forces. In particular, we can mention [START_REF] Bishop | The use of the slip circle in the stability analysis of slopes[END_REF], who supposes that the sliding line is circular and that the vertical component of the inter-slice forces applied to a slice cancels itself out.

CRITERIA OF THE TRAPEZOID FLAT SLIDE

Soft soils depositing within barrage reservoirs are characterized by their low strength and their high water content. These characteristics imply the development of specific slip criteria and the justification of used assumptions. The observation campaigns carried out during several emptying [START_REF] Marot | Stabilité de talus de retenue de barrage[END_REF] enabled us to underline the presence of flat sliding that follows the level of water along the banks of the river. We did not observe discontinuity in the transverse direction. The two dimensional model of flat sliding therefore appears to be adequate to describe in-situ slides.

Thus, we develop a flat-sliding criteria for a slice of soil with a trapezoidal shape. Fig. 1 represents the slice studied, with the forces applied to it.

High level

θ i H H i+1 i P wi Low level h i+1 f i+1 Φ i v i h i v i+1 W i f i Τ ti Σ ti α i b

Fig. 1 Schematic representation of the trapezoid flat slide

The notations used are the following: H i height of water above the side surface i h i height of the side surface separating slice i from the slice i-1 W i weight of slice i P wi force resulting from the pressure of the water applied to the upper surface of slice i Τ ti and Σ ti tangential and normal components of the force applied to the lower surface of the slice i Φ i force applied by slice i-1 to slice i f i and v i horizontal and vertical components of the inter-slice force Φ i α i inclination of the lower surface of slice i θ i inclination of the upper surface of slice i b width of the slice Assumptions: the material is assumed to be a1) homogeneous, a2) saturated, a3) normally consolidated (overpressure is zero compared to the surrounding pressure).

Assumption b): the failure occurs in undrained conditions. Justification of assumption a1: the profiles of total stress and undrained shear strength obtained by different researchers [START_REF] Mekerta | Etude des propriétés géomécaniques des sédiments d'envasement de la retenue du barrage de Génissiat[END_REF][START_REF] Alexis | Mud bank sliding in a dam reservoir[END_REF] for this type of soils allow us to suppose the sediment massif to be homogeneous.

Justification of assumption a2: given the low incidence of non-saturation on the effective stress [START_REF] Marot | Sedimentological and geotechnical approach of correlations between sediment parameters in a reservoir[END_REF] we can assume that the sediment mass is saturated.

Justification of assumption a3: sedimentation and sediment consolidation background storage reservoir is often regular. We can thus suppose that the sediments are normally consolidated prior to emptying. The consolidation phenomenon is often caused by a loading, which creates a variation in the volume of the soil. There is thus a change from a state of hydrostatic balance to a secondary state of hydrostatic balance. For this to happen, there is a dissipation of pore pressure that is accompanied by an increase in the effective stress.

On the other hand we have seen that the mud can be considered as a saturated fluid, thus the lowering of the water level is instantaneously transmitted to pore pressure (incompressible fluid) which remains hydrostatic. Consequently the lowering does not bring about a variation of volume within the material, and does not create pore overpressure. The effective stress thus does not undergo any modification and the material remains normally consolidated. Justification of assumption b: to carry out laboratory drained shear tests on fine soils samples of low thickness, the shear velocity must be very low (lower than 25 µm per min or 3,6 cm per day, AFNOR NF P 94-071-2). However in-situ shearing involves great thicknesses of sediments, and the speed of the drawdown of the storage reservoir, which heavily affects the shear velocity, is generally about several tens of centimetres per day (with a maximum of 10 cm per hour). We can therefore suppose that shearing appearing at the time of the slope slide occurs in undrained conditions.

Expressing Sliding Criteria as Equations

For convenience and to simplify, we consider slice number 1. In accordance with the assumptions of [START_REF] Bishop | The use of the slip circle in the stability analysis of slopes[END_REF], we suppose that the vertical components of the inter-slice forces are equal: ∆v = v1 -v2 = 0. In a sediment mass, the effective horizontal (σ' h ) and vertical (σ' v ) stresses are linked by the earth pressure coefficient (K 0 ):

σ' h = K 0 σ' v (2)
The interstitial fluid is considered saturated, its pressure is noted as u w .

Therefore

σ h = K 0 σ' v + u w γ γ - 1 γ γ K σ σ w w 0 v h + ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =
(3) γ w : density of water (9.81 kN/m 3 ). Whatever the value of K 0 , the total stress ratio in case of soft soil is around the unit. The values of K 0 for the mud being very difficult to measure but around the unit [START_REF] Schlosser | Eléments de mécanique des sols[END_REF], we can suppose that, as in a fluid, σ h = σ v .

The force f 1 can then be expressed by:

L dz z γ h H γ f 1 1 w 0 1 h1 + = ∫ (4) i.e.: h H γ h f 1 1 w 2 1 1 L γ 2 1 + = (5) thus ( ) ( H h H h γ h h 2 2 1 1 w 2 2 2 1 - L - L γ 2 1 ∆f + = ) (6)
With the assessments of the vertical and horizontal forces, we obtain the expression of the safety factor :

( ) θ - α sin α cos f - α sin W α cos b L F P c w u + ∆ = (7)
The study of the stability of the slice is carried out according to two states of loading:

-with the slice completely submerged (high water level) -with the slice completely emergent (low water level)

Sediment massif completely submerged:

We combine Equation ( 6) and:

θ cos L b 2 H H γ P 2 1 w w ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + = (8) ( ) L b γ 2 1 W h h 2 1 + = (9)
to obtain:

( ) ( ) ( ) [ ] cosα 2 α cos θ tan α sin b α cos b 2 F H h H h H H γ D c 2 2 1 1 2 1 w 1 u - - - + + = (10) with, ( ) ( ) [
]

α cos α sin b γ h h h h D 2 2 2 1 2 1 1 - - + = ( 11
)
Sediment massif completely emergent:

In this case, ( )

h h 2 2 2 1 L γ 2 1 ∆f - = (12)
and P w = 0 yield to:

D c 1 u α cos b 2 F = (13)
with D 1 defined by Equation ( 11). Both approaches, mud flow and flat slide with rectangular slice in the undrained conditions [START_REF] Alexis | Mud bank sliding in a dam reservoir[END_REF], thus lead to the same expression of the sliding condition. The mud flow model has the advantage of giving the profile of the movement rates of the fluid, but the water-sediment interface and the shear plane are supposed parallel. The trapezoid criteria is more in conformity with the different configurations that can exist in nature.

Method of Criteria Validation

With the aim to check the validity of the sliding criteria, we have to compare its results with the changes in the thickness of sediment layer during emptying.

This comparison may be carried out in three steps:

-by considering the sediments displaced during emptying, -by considering the sediments not displaced, -by using the criterion within the framework of a displaced height forecast.

The three kinds of surveys are thus referenced (level N, inclination α):

-top level of soft soil before emptying: N 0 , α 0 -top level of soft soil after emptying : N 1 , α 1 -base of soft soil : N T , α T -height of the displaced layer: h d = N 0 -N 1 -height of the soft soil layer: h T = N 0 -N T Fig. 2 summarizes these notations.

To facilitate calculations, we geometrically represented the trapezoidal slice studied (Fig. 3). We propose an overall (spatiotemporal) indexing (n), corresponding to a configuration that may be (according to the case in question), temporal (instant t) or spatial (top level of soft soil). Geometrically, we obtained the following expressions:

h 1 = Z 1 n -T 1 h 2 = Z 2 n -T 2 N 0 N 1 h d α 0 h T α 1 N T α T Fig. 2 Summary of the notations Z e n H 2 n H 1 n Z 2 0 Z 1 h 0 Z 2 n 2 θ Z 1 n h 1 T 2 T 1 α b Fig. 3 Geometrical representation of the trapezoid ( ) T T b 1 2 2 2 b α cos - + = ( ) T T b T T 1 2 2 2 1 2 - α sin - + = ( ) Z Z b n 1 n 2 2 2 b θ cos - + = ( ) Z Z b Z Z n 1 n 2 2 2 n 1 n 2 θ sin - + - =
The heights of water above the top of the trapezoid at the n instant are:

H 1 n = Z e n -Z 1 n H 2 n = Z e n -Z 2 n
We should also emphasize that shear strength in the undrained condition along the potential line of sliding defined by the elevations T 1 and T 2 , can be expressed by: ( )

2 a T Z T Z c c 1 0 1 2 0 2 u0 u - + - + = (14)
-Submerged sediment massif:

Using the preceding expressions and Equation ( 10), we can express the safety factor of the fully submerged sediment massif by:

( ) ( ) [ ] D D T T b c T Z T Z F 2 1 2 2 u0 1 0 1 2 0 2 n 1 2 2 a 2 + + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + - + - = - (15) with, ( ) ( )
( ) ( ) [ ] T Z T Z T T T Z T Z D 1 n 1 2 n 2 2 2 1 2 1 n 1 2 n 2 1 b γ - - - + - - + - = (16) ( )( ) ( )( ) ( ) ] Z Z T T Z Z Z T Z Z Z T Z Z Z [ γ D 2 b 2 n 1 n 2 1 2 n 1 n 2 n e 1 n 1 n 1 n e 2 n 2 n 2 n e w 2 + - - ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - + - - - - - = (17)
-Emergent sediment massif:

Equation ( 13) enables us to express the safety factor in the form:

( ) ( ) [ ] D T T b c T Z T Z F 1 2 2 u0 1 0 1 2 0 2 n 1 2 2 a 2 - + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + - + - = (18)
D 1 is given by Equation ( 16).

(a) Calculation of the safety factor of the displaced sediment

For this study, we need to characterize the sediments located between the level N 0 (top level of soft soil before emptying) and the level N 1 (top level of soft soil after emptying) (Fig. 2). Consequently we assumed that in configuration n, no sliding has occurred yet, i.e. that:

Z 1 n = Z 1 0 and Z 2 n = Z 2 0
These heights correspond to the level N 0 . The potential line of sliding is given by the level N 1 , the heights T 1 and T 2 thus correspond to this level. Before emptying, the safety factor is calculated using Equation (15) (submerged sediment mass) and the preceding conditions. The values for the safety factor thus obtained enabled us to check the stability of the submerged sediment mass.

After emptying, the same conditions were applied to the Equation ( 18). We obtained the safety factor noted F A , making it possible to identify the unstable zones of the exposed sediment mass (F A < 1).

(b) Calculation of the safety factor of the sediments remaining after emptying

The soil remaining is likely to slide on the surface corresponding to the base of soft soil, and its watersediment interface is provided by post-emptying level (N 1 ). The heights T 1 and T 2 are thus represented by level N T , the heights Z 1 n = Z 1 1 and Z 2 n = Z 2 1 correspond to the level N 1 , and the initial top soft soil level N 0 is used to obtain the values of the heights Z 1 0 and Z 2 0 . The safety factor thus obtained is referred to as F B and the values that are lower than the unit make it possible to locate the unstable zones of the remaining sediments.

(c) Calculation of forecast minimal displaced heights

For the unstable slices, the remaining layer must be stable and thus conform to the expression of the safety factor of the emergent sediment mass (Eq. 18), we obtain, then, the quadratic equation:

( ) ( ) ( ) ( ) [ ]0 2 a b γ 2 T T b c T Z T Z F T T Z Z T T Z Z 1 2 n 1 2 2 u0 1 0 1 2 0 2 n 1 2 n 2 n 2 1 2 n 1 2 = + ⎥ ⎦ ⎤ ⎢ ⎣ ⎡ + - + - + - - - + - - (19)
that we can solve, by adopting only the physically acceptable solution:

( ) ( ) ( ) ( ) ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎣ ⎡ + - - - + - + = - + T T b c F T T Z Z T T T T Z 1 2 1 2 2 2 u n 1 2 n 2 n 2 2 1 2 1 n b γ 8 4 2 1 (20) with, 2 a T Z T Z c c 1 0 1 2 0 2 u0 u - + - + =
The heights T 1 and T 2 are given to us by the base level of soft soil and the heights Z 1 0 and Z 2 0 correspond to the top level of soft soil before emptying (N 0 ). The safety factor at limit state is set to the unit. Equation ( 20) thus makes it possible to calculate the elevation of the low point Z 1 n , knowing the elevation of the high point Z 2 n . In order to forecast the final profile, this Z 2 n elevation is unknown and depends on the evolution of the whole profile. To break free from the interaction of the slices, we can make a local assumption on Z 2 n . We chose to consider a particular configuration where, for a given slice, the slide results in the lowering of the low point, the high point remaining fixed. The minimal displaced height (h dmin ) taken into account is the lowering of the mid point of the slice, that is to say:

Z 2 n = Z 2 0 and ( ) Z Z h n 1 0 1 dmin 2 1 - = (21)

CONFRONTING THE CRITERIA WITH BATHYMETRIC CHANGES

Description of the Selected Site

The storage reservoir studied is that of the Grangent dam which is located on the Loire (France) close to Saint-Étienne. Built between 1955 and 1957, the structure is a cylindrical concrete vault with a height of 55 m and a crown length of 200 m.

The storage reservoir has a volume of 57.4 Mm 3 , a length of 23 km and the surface area of the total catchment is 3850 km 2 . Into this storage reservoir, the Loire receives three tributaries, these being the Semène, the Ondaine and the Lizeron. The study zone extends upstream and downstream from the confluence of the Ondaine with the Loire, for approximately 4 km. We distinguish four zones: upstream, Ondaine, post-confluence and downstream (Fig. 4).

Fig. 4 Positioning of the sampling and bathymetric profile stations

Topographic Profiles

Fifty land surveys were carried out every 400 m in 1952, i.e. before the construction of the dam. These make it possible to know the elevation of the base of the soft soil (N T ).

Bathymetric Profiles

During the emptying between November 1995 and May 1996, the storage reservoir was brought back to a height of 404 m NGF, this being 16 m below its nominal elevation. This lowering of the water table reduced by half the length of the storage reservoir, whose tail then reached the confluence with the Ondaine. An objective being to highlight the influence of the lowering of the water level on the geomorphology of the upstream slopes, two bathymetric survey campaigns using echo sounder were carried out by the Laboratoire National d'Hydraulique (EDF). The first batch of surveys were taken before partial emptying (level N 0 ) and the second batch upon completion (level N 1 ). Positioning was achieved by satellite (GPSD) in order to improve the accuracy of a comparative study of the profiles. Over half of the surface area of the sediments underwent displacement due to the flow of the river or slide on the banks. About the quarter of the sediments present in the storage reservoir moved in this way, to settle in the post-confluence area. The thickness variation of the silt deposition is very marked :

-upstream area (22 profiles): from -0.24 m to -2.40 m and -1m on average, -post-confluence area and downstream (10 profiles): up to +3.40 m and +1.70 m on average.

-

The weakness of the suspended solids concentrations in water downstream of the dam shows that this deposition occurs totally before reaching the structure. Most of the sediment is then concentrated in the post-confluence area and downstream.

Samples and Associated Measurements

With the aim of characterizing all of the sediments submerged before emptying, a sampling campaign associated with laboratory measurements is carried out. Fourteen core samples of a diameter of 6 cm and an average length of 2.5 m were taken equitably over the whole surface of the pre-defined study area to detect possible heterogeneities of the sediments. Several samples were taken at each sampling point in order to eliminate singularities non-representative of the true state of the soil.

We carry out a multivariate analysis in order to define all variables according to three factors and to thus underline the presence of correlation relationships between the various values [START_REF] Marot | Sedimentological and geotechnical approach of correlations between sediment parameters in a reservoir[END_REF]. The correlation analysis shows that the profile representative of undrained shear strength may be estimated by: c u = c u0 + a h c u0 = 0.16 kPa a = 1.47 kPa/m (r = 0.858; n = 48; ξ = 8.9 10 -11 )

with r: correlation coefficient; n: number of individuals; ξ: risk of error

The calculation of the total (σ) and effective stresses (σ') according to the vertical elevation in the sediment (h) leads to:

σ [kPa] = 10.17 h [m] (r = 0.994, n = 177, ξ < 1.10 -11 ) (23) σ'[kPa] = 2.71 h [m] (r = 0.927, n = 177, ξ = 1.10 -11 ) (24)

Locating the Unstable Zones and Calculation of Minimal Displaced Height

The eleven profiles for which all of the bathymetric and topographic data are accessible, and at the level of which the mud mass became exposed during emptying are studied. The values of the safety factor of the submerged sediments (Eq. 15) being all higher than the unit, the model confirms the stability of the sediment masses submerged before emptying. The values lower than the unit of the safety factor (F A ) of the displaced sediments indicate the location of the unstable zones created during emptying.

Fig. 5 shows these unstable zones for one of the profiles studied (P32). These zones, which characterize the sediments before emptying, are thus located on the slopes of the profile N 0 . We note that the central part of the profile, only slightly sloping, is stable with the slide but undergoes major erosion, (not taken into account in this criteria). This same remark can also be made for the profiles located far upstream in the storage reservoir, for the slopes are very gentle there. On level N 1 , we can see the zones of instability of the remaining sediments, characterized by a value for the safety factor F B lower than the unit. We can thus note that sliding has made it possible for the massif to recover its stability. On this graph, all of the points with a x-value equal to zero correspond to the central parts of the profiles, which are the places of slight sloping and of major erosion. The other points, corresponding to the embankment, show a calculated minimal height slightly lower than the measured height. This graph thus presents a satisfactory estimate of the heights displaced by sliding. The trapezoidal shape allows to obtain better estimations than with a rectangular shaped slice.

CONCLUSION

The aim of this study is to simulate mud bank sliding during a dam reservoir emptying process. A flatsliding geotechnical criteria is developed for a slice of soil with a trapezoidal shape. This sliding criteria is based on the theoretical formulation of the safety factor when the massif is submerged and emergent. The data necessary to the criteria are available before emptying.

This criteria first of all made it possible to check the stability of the submerged sediment masses and to locate the unstable zones which appear as the mud is exposed. Then we check that the slide allow for the remaining sediments to recover their stability.

Lastly, with prediction in mind, we estimate the mud thicknesses which would be displaced during emersion in order for the remaining sediments to remain stable.

These estimates are shown to agree with the measurements carried out on the selected site. In the zones subjected to heavy erosion, the thicknesses displaced are underestimated.

However, the simplicity of this criterion means that it can be coupled with a model of erosiondeposition in order to better allow for the representation of the true mechanisms of the soft soil movement caused by the emptying.

  the storage reservoir in configuration n, T 1 and T 2 : elevations of the shear plane.

Fig. 5

 5 Fig.5Localization of the zones of instability, before and after emptying (profile P32) Fig.6makes it possible to compare all the measured displaced heights (h d ) and minimal heights calculated by the criteria, for all of the studied profiles.

Fig. 6

 6 Fig. 6 Measured displaced height vs calculated minimum displaced height

Doctor,

Doctor-Engineer,

Professor, IUT de Saint Nazaire, Laboratoire Génie Civil Nantes-Saint Nazaire, BP

420, 44606 SAINT NAZAIRE cedex, FRANCE, Tel.: 00 33 2 40 17 81 50, Fax: 00 33 2 40 17 81 60,