Maxime Romanet
email: maxime.romanet@femto-st.fr

Luis Miguel

Giraldo Salazar

Maxime Zerbib

Etienne Rochat

Kien Phan Huy

Jean-Charles Beugnot

Luis Miguel Giraldo

Towards single-photon Brillouin optical time domain reflectometry

Introduction

Brillouin scattering based distributed optical fiber sensors have become essential tools for civil engineering applications such as monitoring of pipelines or submarine cables [START_REF] Motil | State of the art of brillouin fiber-optic distributed sensing[END_REF][START_REF] Hartog | An Introduction to Distributed Optical Fibre Sensors[END_REF]. Their field of application has expanded opening up new horizons for geophysicists and natural scientists [START_REF] Wang | Test on application of distributed fiber optic sensing technique into soil slope monitoring[END_REF][START_REF] Yin | Real-time monitoring and early warning of landslides at relocated Wushan Town, the Three Gorges Reservoir, China[END_REF][START_REF] Zeni | Brillouin optical time-domain analysis for geotechnical monitoring[END_REF]. These sensors, which take advantage of the frequency dependence of Brillouin gain on temperature and strain, fall into two categories [START_REF] Hartog | An Introduction to Distributed Optical Fibre Sensors[END_REF][START_REF] Bai | Recent Advances in Brillouin Optical Time Domain Reflectometry[END_REF]. Brillouin optical time domain analysis sensors (BOTDA) use stimulated Brillouin scattering and require access to both ends of the fiber. Brillouin optical time domain reflectometers (BOTDR) that measure spontaneous Brillouin scattering require access to only one end of the fiber and provide a clear advantage in terms of simplicity. Increasing signal-to-noise ratio (SNR) represents a major challenge for these sensors. It determines the range of the sensor, typically 75 km on commercial devices. Increasing the range of BOTDR sensors makes it possible to contribute to the ecological transition, by helping to guarantee the integrity of submarine electrical cables for offshore wind farms. Numerous improvements have been proposed in terms of probe light [START_REF] Li | SNR enhancement in Brillouin optical time domain reflectometer using multiwavelength coherent detection[END_REF][START_REF] Zhang | Performance improvement for long-range BOTDR sensing system based on high extinction ratio modulator[END_REF], remote amplification with Raman [START_REF] Cho | 50-km single-ended spontaneous-Brillouin-based distributed-temperature sensor exploiting pulsed Raman amplification[END_REF][START_REF] Cho | Enhanced performance of long range Brillouin intensity based temperature sensors using remote Raman amplification[END_REF][START_REF] Song | 100 km Brillouin optical time-domain reflectometer based on unidirectionally pumped Raman amplification[END_REF] or EDFA amplifiers [START_REF] Gyger | Extending the sensing range of Brillouin optical timedomain analysis up to 325 km combining four optical repeaters[END_REF][START_REF] Clement | Enhancement of sensing range of Brillouin optical time-domain reflectometry system up to 150 km with in-line bi-directional erbium-doped fibre amplifications[END_REF] or coding strategies [START_REF] Soto | Analysis of optical pulse coding in spontaneous Brillouin-based distributed temperature sensors[END_REF][START_REF] Soto | Brillouin-Based Distributed Temperature Sensor Employing Pulse Coding[END_REF][START_REF] Wang | Enhancing the performance of BOTDR based on the combination of FFT technique and complementary coding[END_REF]. The current range of BOTDR sensors is 100 km with Polarization-Diverse Coherent Detection [START_REF] Jostmeier | Long-distance botdr interrogator with polarization-diverse coherent detection and power evaluation[END_REF], and 250 km combining Raman amplification and an EDF after few kilometers [START_REF] Joy | Increasing the distance range of repeaterless Brillouin -otdr to 250 km by optical amplification[END_REF]. The former achieves 100 km temperature measurement with a straightforward design but also a signal-to-noise ratio degraded after 80 km. The latter, along with most recent improvements involves adding amplifiers within the fiber network or added complexity to the BOTDR setup.

Eventually, the signal-to-noise ratio (SNR) of the photodetector always limit the ultimate performance. With the development of quantum technologies, new detectors such as single-photon avalanche diodes (SPADs) and superconducting nanowire single-photon detectors (SNSPDs) become available [START_REF] Eraerds | Photon Counting OTDR: Advantages and Limitations[END_REF]. They have shown their efficiency in many OTDR applications, improving sensing range while maintaining a good resolution [START_REF] Diamanti | 1.5 µm photon-counting optical time-domain reflectometry with a single-photon detector based on upconversion in a periodically poled lithium niobate waveguide[END_REF][START_REF] Hu | Photon-Counting Optical Time-Domain Reflectometry Using a Superconducting Nanowire Single-Photon Detector[END_REF][START_REF] Zhao | Long-haul and highresolution optical time domain reflectometry using superconducting nanowire single-photon detectors[END_REF][START_REF] Shentu | 217 km long distance photon-counting optical time-domain reflectometry based on ultra-low noise up-conversion single photon detector[END_REF]. They have also been used in distributed temperature sensors (DTS) using Raman backscattering, with a range of 50 km [START_REF] Antony | Single-photon detector based long-range fibre-optic distributed temperature sensor[END_REF], and achieve spatial resolution of the order of a centimeter [START_REF] Tanner | High-resolution single-mode fiber-optic distributed raman sensor for absolute temperature measurement using superconducting nanowire single-photon detectors[END_REF][START_REF] Gasser | Distributed temperature sensor combining centimeter resolution with hundreds of meters sensing range[END_REF]. They are thus good candidate to help developping a class of new devices, single-photon Brillouin optical time domain reflectometers (ν-BOTDR), that could complement the traditional BOTDRs.

However, few studies have been carried out on ν-BOTDR. They use either the Landau-Placzek ratio (LPR) requiring Rayleigh and Brillouin intensity measurements, to determine the temperature [START_REF] Xia | Brillouin optical time-domain reflectometry using up-conversion single-photon detector[END_REF] or a scanning Fabry-Perot interferometer to scan the Brillouin resonance [START_REF] Xia | A Distributed Brillouin Temperature Sensor Using a Single-Photon Detector[END_REF]. Measurements were reported up to 50 km with an estimated range of 90 km [START_REF] Xia | A Distributed Brillouin Temperature Sensor Using a Single-Photon Detector[END_REF]. Therefore, numerous questions are left to be addressed to build a ν-BOTDR such as range improvement, retrieving temperature information or ability to discriminate temperature from strain.

In this paper, we investigate a single-photon ν-BOTDR setup to perform distributed temperature measurements using a narrow fiber Bragg grating (FBG) without frequency scan. This idea has already been used in schemes using traditional detectors with Mach-Zehnder interferometers as discriminators [START_REF] Masoudi | Distributed dynamic large strain optical fiber sensor based on the detection of spontaneous Brillouin scattering[END_REF]. In particular, our configuration allows any instability issue to be addressed with a simple post-processing of the data. In the first section, we will present our experimental setup and demonstrate its detection range. Then, we will explain the principle of temperature measurement with a narrow FBG used as a frequency discriminator and show our experimental results. In the third section, we will show how instabilities can be corrected through real-time monitoring and post-processing of the data, thus eliminating the need for a feedback loop. Finally, we will give an overview of how our setup can deal with temperature-strain cross-sensitivity.

Measurement setup and sensing range

Our experimental setup is shown in Fig. 1 (top). A tunable narrowband laser (Apex technologies) sends light signal at λ = 1547 nm into a polarization controller (PC). The light is split in two beams by a 90/10 coupler. The beam with 10% of the light reaches an optical isolator (ISO) and a semiconductor optical amplifier (SOA) that shapes the pump pulses. Pulses go through a second optical isolator and an optical circulator (C1) before entering the fiber under test (FUT): a 120 km single-mode fiber (G625D). In order to carry out measurements automatically, all components are controlled by a computer. Within the FUT, the light experiences spontaneous Brillouin scattering. It is scattered by thermally excited phonons that propagate along the fiber. The back-scattered light is shifted in frequency by the Doppler effect. The frequency shift is related to the acoustic velocity in the optical fiber, which depends on the temperature and the strain. The back-scattered light returns to the circulator C1 where it is routed towards an optical switch SW1. In measurement mode, the backscattered light is sent through SW1 to a circulator C2 and a narrow Fiber-Bragg Grating (FBG) that serves as a frequency discriminating reflector. Then the filtered backscattered light reaches a second switch SW2. SW2 carries the backscattered photons to an infrared SPAD from Auréa Technology operated in the so-called gated mode or Geiger mode. The applied bias is above the breakdown voltage, in order to have sufficient gain to be able to detect a single incident photon. The SPAD is based on InGaAs/InP, cooled around -50°C.

In order to make a distributed measurement along the optical fiber, it is important to locate precisely the position at which the light was backscattered; this is achieved by a time-of-flight measurement. As shown in Fig. 1 (top), by controlling the time delay ∆t between the launch of the pulse from the SOA (blue 1) and the activation gate of the SPAD (pink 2), we can carry out measurements at a particular distance d:

∆ t = 2d c n SMF + t setup ( 1
)
where c is the speed of light in vacuum, n SMF is the fiber refractive index, d is the sensing distance, t setup is the delay due to propagation time in setup components (isolator, circulators, switch, patch cords, filter) and ∆t is the photon round trip time of flight in the fiber under test (FUT).

The main advantage of using SPADs over conventional detectors is to increase the sensing range of the distributed sensor, which is determined by the signal-to-noise ratio (SNR) of the whole setup. In our case, the SNR is the ratio between the measured Brillouin anti-Stokes signal and the noise level on the SPAD. The strength of the Brillouin anti-Stokes signal is determined by the peak power of the SOA pulses, the losses in the optical fiber and components, the Brillouin scattering efficiency and the detection efficiency of the SPAD. The noise on the SPAD is made up of two types of noise, the optical noise from the setup, and the intrinsic noise of the SPAD. The optical noise is related to the optical signal coming from the setup. It is determined by three factors. First, the extinction rate of the SOA defines the background noise between the pulses. This noise can eventually reach the SPAD and trigger unwanted photon counts. Then, the cross-talk of the circulators C1 and C2 quantifies the signal that leaks from the SOA directly to the SPAD without going through the fiber under test. Finally, the filter efficiency gives the rejection efficiency of the Rayleigh, Brillouin Stokes and Raman backscattering. The second is the intrinsic noise of the SPAD, called dark count rate (DCR), which depends on the efficiency, the dead time and the duration of the gates set on the SPAD. As a consequence, careful attention has been given to the setup to provide the best SNR with the available components, and also to the SPAD parameters chosen in order to limit the DCR and avoid after pulsing. These parameters are specified in the table in Fig. 1 (bottom).

Using the specifications in Fig. 1 (bottom), we can predict the detection rate:

C ,cps = η s τ T S B,W ω , (2)
Where is the reduced Planck constant, ω is the pulsation of light, η s the SPAD efficiency, τ the gate duration of the SPAD, T the period and S B,mW the backscattered Brillouin power [mW] defined as: The backscattered Brillouin power along the fiber is defined as:

S B,mW = 10 S B,dBm
S B,dBm = P peak -g B -2 d α SMF -α setup , [dBm] (4)
where P Peak is the pulse peak power injected in the fiber, g B is the Brillouin efficiency (78 dB in our configuration), d is the probed distance [km], α SMF are fiber losses and α setup are components losses of the setup (circulators, filter, switch).

The experimental results are presented in Fig. 2. The Brillouin anti-Stokes count rate is reported in count per second (cps) with a logarithmic scale as a function of the probed distance.

The low noise level of the sensor (0.05 cps), represented in blue on Fig. 2 allows us to reach a distance up to 120 km, with a 10 m spatial resolution, which is 30 km longer than the results reported in [START_REF] Xia | A Distributed Brillouin Temperature Sensor Using a Single-Photon Detector[END_REF]. On Fig. 2 we only measure from 50 km where current commercial devices start to be less efficient. These results could be further improved by increasing the signal-to-noise ratio (SNR) of the setup. For example, only +15 dBm power is launched into the fiber while commercial devices commonly use +23 dBm. Similarly, using a SNSPD, one could improve the SNR by 20 dB [START_REF] Quantique | ID281 Superconducting Nanowire Series[END_REF]. In theory with proper components such a setup could then reached over 190 km, approaching the performance obtained with a ν -OT DR, using a SNSPD, with a range of 217km [START_REF] Shentu | 217 km long distance photon-counting optical time-domain reflectometry based on ultra-low noise up-conversion single photon detector[END_REF]. The difference between the theoretical range of our setup and the range of ν-OTDR is simply due to the higher efficiency of Rayleigh backscattering. In this article, at 120 km, the sensor still has a signal-to-noise ratio (SNR) of 3.6 dB, which is a first step in the demonstration of a ν-BOTDR.

Temperature measurements

BOTDR sensors use the Brillouin frequency shift, to measure the temperature along an optical fiber. Usually, a heterodyne setup with a local oscillator is used. Such a configuration is not allowed in our case because the local oscillator would blind our single-photon detector. In our experiment, we develop a frequency discriminator. This allows us to transform a frequency variation into an intensity variation. The frequency discriminator is the slope of a narrow FBG that gives access to the Brillouin frequency. The FBG frequency response is shown in Fig. 3, where the linear area of the slope (over 576 MHz) is represented in red. From the frequency shift, we can deduce the temperature shift with the temperature sensitivity coefficient C ν B T = 1.07 MHz/K for a step index single-mode fiber at 1550 nm [START_REF] Hartog | An Introduction to Distributed Optical Fibre Sensors[END_REF]. Thanks to this method, we only need to set the frequency laser to the desired working point shown in Fig. 3, and during the measurement we do not need to make a frequency scan of the Brillouin response. By avoiding the frequency sweep, we also reduce the measurement time.

It shows that, for an operating point halfway down the left-hand slope of the transmission band, we get 5% transmission variation for 42 MHz frequency shift. This ratio η c = 8.2 MHz/% enables to turn a frequency shift variation into a count rate variation on the SPAD. In order to perform distributed temperature measurement, the propagation losses shown in Fig. 2 must be removed. Count rate datas are thus processed in order to compensate for the fiber losses, as a function of the probed distance, as depicted in Fig. 5.

In Fig. 4 (a), we report the temperature measurement obtained this way at operating point A (see Fig. 3) from a 10 m sample of a fiber immersed in 85 ± 1.5 • C water after 100 km.

We demonstrate the ability to detect a hot point of 10 m, at 100 km, with 100 ns optical pulses. The profile is realized with a sampling interval of 2.5 m. In gated mode, the sampling interval corresponds to the time delay between two active gates on SPAD. At this distance, we get a count rate of 0.55 count per second (cps), cor- responding to a SNR of 10.4 dB. The points measured at 45 °C outside the bath simply correspond to samples overlapping the thermal transition, linked to a spatial resolution of 10 m. Moreover, the jitter delay line of 1.5 ns adds uncertainty to the sampling. In Fig. 4 (b), we present the standard deviation on 10 consecutive temperature measurements according to the acquisition time at 60 km and 100 km. We demonstrate experimentally a standard deviation of σ 100km = 8.9 • C at 100 km and σ 60km = 5.4 • C at 60 km, with an acquisition time of 10 minutes for both distances. The difference of the standard deviation between the 2 distances is explained by the higher level of signal at 60 km (+15.2 dB compared to 100 km), giving a better statistic and higher number of events. One can also see that for sufficient time acquisition, the standard deviation reaches an optimum that defines the temperature resolution. Even if we use the FBG slope to avoid a frequency scan, the acquisition time is quite long compared to conventional BOTDR. This is due to an acquisition point per point along the optical fiber which is intrinsic to the gated mode of the SPAD. This measurement confirms the narrow FBG as a valid frequency discriminator, to convert count rate variation into temperature measurements. A more detailed reduction of acquisition time measurement of our ν-BOTDR is the subject of a future work.

Real-time monitoring and data post processing

The use of a frequency discriminator such as a narrow FBG has an obvious disadvantage. The relative frequency deviation between the laser and the FBG frequency is expected to be perfectly stable. Such setup will be sensitive to any form of instability (laser drift, high temperature variation). It thus usually requires a feedback loop to be reliable at the cost of increased complexity. However, this is not necessary in our case if we take advantage of the statistical nature of single-photon measurements. Indeed, the conversion of the count rate into Brillouin frequency shift is done by knowing the slope of the FBG at the given operating point. This operating point varies with the relative frequency difference between the laser and the FBG frequency. By measuring this difference, a new operating point can be deduced and the Brillouin frequency shift calculated with an updated conversion parameter η c . Since single-photon detectors record events, it is only necessary to apply the measured relative frequency position of the FBG on experimental data in order to correct the drift as shown in Fig. 5. To achieve this, first we use a FBG with an athermal package in order to prevent thermal shocks and rapid temperature variations, leaving only slow temperature drifts (>10 min). Every 10 minutes, the setup switches to monitoring mode and the relative deviation between the laser and the FBG is monitored using the part of the setup highlighted in green in Fig. 1. Switches SW1 and SW2 are then used to take some of the light from the laser and pass it through an EOM modulator. The modulator produces a frequency-adjustable harmonic that probes the FBG filter at 11 GHz.

A frequency sweeping of the laser allows scanning the FBG, and the transmission is measured by the photodiode PD. Thus, FBG frequency shift with respect to the laser frequency is measured. The position of the operating point at 11 GHz of the laser frequency for a temperature of 300 K can be computed (red area on Fig. 3). If the FBG has drifted with respect to the laser, the slope at the operating point is measured and the following registered single-photon detections will be post processed accordingly as shown in green in Fig. 5. Thanks to the FBG scanning, it is possible to compensate for filter drift, and thus to correct frequency measurements thanks to a dedicated post-processing. Note that this does not add complexity to the detection setup since single-photon measurement is inherently statistical. Single-photon detection requires the electronics to record an event for each detected photon and then calculate the average photon flux. Adding additional processing to this calculation step is of negligible cost, while the lack of a feedback loop makes the setup simple and robust in measurement mode. Standard BOTDR setup usually requires an optoelectronic setup to continuously scan the Brillouin frequency. Here the optoelectronic setup is only used to monitor the system. Similarly, any other environment monitoring measurement could be added to retain or discard specific recorded events and refine the measurement.

Towards cross-sensitivity measurements

One of the caveats of our setup is that the count rate we measure on working point A (see Fig. 3) depends on both the Brillouin frequency shift and the efficiency of the Brillouin scattering. Both varies with temperature and strain [START_REF] Hartog | An Introduction to Distributed Optical Fibre Sensors[END_REF]. One way to overcome this problem is to perform an additional control measurement when a temperature variation is recorded. By setting the laser frequency to change the operating point to B in the middle of the FBG band pass filter as shown in Fig. 3, we only measure the efficiency of the Brillouin scattering. The superposition of both measurement are reported in Fig. 6, where the sensitivity is computed for the two methods. One can see in red the measurement with operating point A that includes both frequency shift and intensity contribution, while the measurement in blue only gives the intensity contribution. On Fig. 6 we show the sensitivity of both methods, defined as:

Sensitivity = Count rate avg(Count rate Ambiant ) × 100 [%], (5)
where Count rate Ambiant corresponds to the count rate averaging between 975 and 980 m, which is the reference level signal at room temperature. We clearly see that sensitivity at working point A, by using the slope of the FBG is higher (41%) than the B working point (16%) even if the count rate is lower for operational point A due to FBG losses. The highest sensitivity is obtained thanks to the frequency sensitivity ratio η c linked to the slope. Thereby, our method demonstrates a better sensitivity compared to LPR method based on Brillouin and Rayleigh intensity measurements [START_REF] Xia | A Distributed Brillouin Temperature Sensor Using a Single-Photon Detector[END_REF].

Similarly to [START_REF] Hartog | An Introduction to Distributed Optical Fibre Sensors[END_REF], if we measure the count rate variation for a given temperature and strain variation for both A and B operating points, we can obtain the coefficients of the matrix M as seen on the right-hand side of:

∆C A ∆C B = a ∆T,A a ∆T,B a ∆ε,A a ∆ε,B ∆T ∆ε , (6)
= M ∆T ∆ε .

By inverting this matrix, one can retrieve the temperature and strain from both count rates variation.

∆T ∆ε = M -1 ∆C A ∆C B (7)
Where ∆C A and ∆C B are respectively the count rate variation measured for working point A and B. The temperature shift coefficients a ∆T,A , a ∆T,B and strain shift variation a ∆ε,A , a ∆ε,B are related to temperature variation ∆T and strain variation ∆ε with [START_REF] Bai | Recent Advances in Brillouin Optical Time Domain Reflectometry[END_REF]. A measurement at the working point A and B would allow performing temperature and strain measurements separately.

Conclusion

We have reported a ν-BOTDR temperature measurement based on SPAD and a FBG frequency discriminator. The advantage of using SPADs for distributed sensor is the ability to measure small signal down to -120 dBm, without the need for a cryostat, unlike SNSPDs. The high signal-to-noise ratio (SNR) allows measuring the distributed Brillouin response up to 120 km with a spatial resolution of 10 m in standard single mode fiber. We demonstrate the ability to detect a hot spot at 100 km, and obtain a higher sensibility using a FBG compared to the LPR method. Thanks to the monitoring mode, environmental instabilities can be taken care of without the need of a feedback loop or adding complexity to the detection setup. By using different operating points on the FBG we show the possibility to discriminate strain and temperature effect without measurement of the Rayleigh intensity. In the future, the signal-to-noise ratio offered by single-photon detectors can open new possibility to extend the measurement range of Brillouin sensor reflectometers. Such performances can be useful for monitoring long submarine power cables, when fiber loop are not available and where no power is available along the fiber to amplify the optical signal. This new device might thus help with the ongoing ecological transition.

Fig. 2 .

 2 Fig. 2. (Red) ν-BOTDR trace from 50 km to 120 km realized with operating point A (see Fig. 3), with a 10 m spatial resolution (SOA 100 ns pulse width), 50 ns SPAD gates. (Blue) Noise floor average. (Dark dotted) Noise floor standard deviation.

Fig. 3 .Fig. 4 .

 34 Fig. 3. Frequency response of our narrow FBG, scanned with SPAD, with an acquisition times equal to 60 s. Red circle: operating point A for frequency measurement. Blue star: operating point B, for intensity measurement only. The portion of the slope highlighted in light red indicates the area in which the response is linear, with a maximum deviation of 10% from linearity. It extends over 576 MHz.

Fig. 5 .

 5 Fig. 5. Bloc diagram explaining the real-time monitoring of the ν-BOTDR. In red, the data processing of the measurement mode is explained. In green, the monitoring mode is shown. It allows scanning the FBG shift every 10 minutes, in less than 1 minute, in order to compute the new operating point (OP) and the FBG slope. The Brillouin anti-Stokes resonance is represented in blue on the diagram. The contribution of the new operating point to the statistical measurement of the Brillouin frequency shift is corrected accordingly.

Fig. 6 .

 6 Fig. 6. Comparison of 2 working points A and B (see Fig. 3), with a 10 m fiber spool in bath at 84 ± 1.5 • C, at 1 km, with 15 ns optical pulses. Blue: Measurements realized at operating point B. Red: Measurements realized at operating point A. (a) Original data measurement. (b) Sensitivity comparison of both methods, with a reference at ambient temperature between 975 m and 980 m by averaging counts on this area.

Acknowledgment

We thank Marc Niklès and Félix Bussières for fruitful discussions, Jean-Marc Merolla for helpful discussion and AUREA technology for support in SPAD.

Funding

The authors would like to acknowledge the financial support of EIPHI Graduate School (contract ANR-17-EURE-0002), European ACTPHAST 4.0 project and French research ministry and SUPMICROTECH-ENSMM.

Disclosures

The authors declare no conflicts of interest.


Guillaume Castera

Juliette Chabassier

Paul Fisette

Brad Wagijo

Piano bridge mobility and longitudinal precursors

Keywords: Piano bridge, Soundboard, Coupling, Numerical simulation, Power balance, String orientation, Transmitted power Chevalet de piano, Table d'harmonie, Couplage, Simulation numérique, Bilan d'énergie, Orientation des cordes, Puissance transmise

This report is devoted to the modelling of the complex movements of the piano bridge for the transmission of waves between the strings and the soundboard. The model allows accurate transmission of longitudinal waves by taking into account the rotational movements of the bridge and their transmission to the shear directions of the soundboard, modelled as a Reissner-Mindlin plate. The addition of these mobilities, which to our knowledge had never been modelled before in the context of the coupling between a non-linear string and a soundboard, makes it possible to observe longitudinal precursors. We also propose a study of the influence of the angular orientation of the string relative to the soundboard on the power transmitted via the bridge. This study, made possible by our numerical simulation tool, could have applications for instrument making.

Introduction

Although the strings are the elements that produce the sound through their vibrations, a piano, like all stringed instruments, would be inaudible without a soundboard, because the surface area of the strings in contact with the ambient air is too small. A soundboard is a thin wooden plate that allows the vibrations of the strings to be transmitted to the ambient air with a much larger surface area and therefore much greater power. The choice of wood for soundboards, its growth, cutting, drying, etc, is part of the expertise of the instrument makers and contributes greatly to the sound produced.

To transmit the vibrations from the strings to the soundboard, an element is needed to create an acoustic bridge between them. This is the purpose of the piano bridge. This is a wooden beam, usually between 20 and 60 millimetres high, glued to the soundboard and over which the strings are attached. Piano makers pay particular attention to the adjustment and construction of the bridge, as it has a considerable influence on the sound [START_REF] Corradi | Modal analysis of a grand piano soundboard at successive manufacturing stages[END_REF]. Through years of evolution, it is one of the most complex mechanisms on the piano, as evidenced, for example, by the zig-zag attachments of the strings to its surface. Although originally made from a single piece of wood, today's bridges are mostly laminated and assembled from wood. This gives greater flexibility in the choice of species, and also reduces the risk of cracking. Some bridges are even drilled to reduce their weight and increase their mobility.

The vast majority of modern pianos have two separate bridges for the bass and treble strings. This is the result of a compromise aimed at maximising the length of the strings while keeping their attachment points as far away as possible from the fixed edges of the soundboard for optimum sound rendering. In fact, just as on a drum, a strike close to the fixed edge emits much less sound than a strike centred on it. To further optimise this point of attachment in relation to the length of the strings, many of today's pianos have a cantilever bass bridge which allows them to be closer to the centre of the soundboard. Finally, the down bearing, i.e. the load exerted by the strings on the soundboard through the bridge, is certainly one of the most important and sensitive adjustments made by piano makers. For optimum sound transmission, neither too much pressure, which would mu e the vibrations, nor too little pressure, which would transmit too few waves, is required.

Soundboard

Cantilever bridge String

Bearing (load)

Down bearing angle

Inria

In practice, this bearing is adjusted by setting the angle of the strings from the horizontal, known as the down bearing angle. At zero angle, the string exerts no pressure on the soundboard and transmits very little vibration. The greater the angle, the greater the force exerted on the bridge and the greater the risk of mu ing the sound [Reblitz, 2005].

The aim of our study is not to reproduce exactly all the complexities present at this point on the instrument, but to model a mechanical behaviour that can account for the transmission of waves, particularly longitudinal waves. Experimental studies have been carried out in [Conklin, 1996] and [Giordano, 1998] and identify two possible explanations for the transmission of longitudinal waves via the bridge: the down bearing angle and/or a rotational and torsional movement of the bridge. Bridge models already exist in the literature and can be used to test some of these hypotheses. For example, [START_REF] Cuenca | Three-dimensional interaction between strings, bridge and soundboard in modern piano's treble range[END_REF] and [Pfeifle, 2014] model the transmission of the transverse wave from a linear Euler-Bernoulli string to a soundboard. [START_REF] Valiente | Modeling the interaction between piano strings and the soundboard[END_REF] focuses on the di erent transversal polarizations of a linear Euler-Bernoulli string coupled to a soundboard through a mobility matrix. The previous models being linear, they do not account for the longitudinal waves. On the contrary, [Tan, 2017] models the interaction of a doubly polarised non-linear string with a bridge that is mobile in all directions in space, but does not allow coupling to a soundboard. [Chabassier, 2012] couples the transverse and longitudinal vibrations of a non-linear string to a bridge movable perpendicular to the soundboard, with a loading angle. The author notes that longitudinal waves are present, but largely underestimated if the angle corresponds to that of a real instrument. To observe correct transmission, the angle must be overestimated and no longer corresponds to a possible configuration on a piano, which implies that another mechanism must allow longitudinal waves to be transmitted. However, there is currently no model that takes into account the rotational or torsional movements of a bridge coupled to a soundboard, or the e ect of these movements on the transmission of longitudinal waves.

In [START_REF] Podlesak | Dispersion of waves in piano strings[END_REF], the authors carry out the following experiment. An accelerometer is placed on the bridge to measure its mobility perpendicular to the soundboard, and a microphone is placed in the room to record the sound. They observe that the microphone detects a sound signal even before the accelerometer detects the presence of a wave at the bridge. This confirms that a longitudinal wave, with a faster propagation speed than the transverse wave, propagates in the string and forms a longitudinal precursor. We would like our model to accurately represent the presence of this precursor.

In this report, we present a model of wave transmission between a string and a soundboard modelled as a Reissner-Mindlin vibrating plate [Reissner, 1945], taking into account the main spatial movements of the bridge. In particular, we add to the model presented in [Chabassier, 2012] the rotational movements of the bridge and their coupling with the shear components of the table. The model and its discrete equivalent preserve power balances to ensure the numerical stability of the coupling. Finally, we show that these refinements allow longitudinal waves to be transmitted and longitudinal precursors to be observed with realistic down bearing angles.

The soundboard as a Reissner-Mindlin plate

Piano soundboards are quite thin, usually 6 to 10 millimeters thick. Consequently, one of the three dimensions (the thickness) is very small compared with the other two (the length and width), which can measure several meters on large concert pianos. This observation led to consider a reduced model of the equations of 3D linear elastodynamics, and to model the soundboard not as a threedimensional solid, but as a plate with only two dimensions. Existing plate models include the Kircho -Love model, adapted for thin plates and which can be compared to the d'Alembert model for a one-dimensional string, and the Reissner-Mindlin model for thicker and sti er plates, whose one-dimensional equivalent would be the Timoshenko model.

The Reissner-Mindlin model

In continuity with the work of [Chabassier, 2012] we use and describe in the following the Reissner-Mindlin model detailed in [Grob, 2006] for which three equations describe a membrane motion of plate elevation u p and two shear motions described by two angles

◊ p = t (◊ 1 p , ◊ 2 p ) on a bounded domain µ R 2 .
" is the plate thickness, fl the density, E x and E y the Young's moduli, G xy G xz and G yz the shear moduli, Ÿ the shear correction coe cient, ‹ the Poisson's ratio and Á the linearized strain tensor. Here, z is the direction perpendicular to the plate. We begin here by writing down the boundary conditions for clamped plate, i.e. the Dirichlet conditions.

Seek (u p , ◊ p ) : ◊ [0, T ] ≠ae R 3 so that (P) Y _ _ _ _ ] _ _ _ _ [ fl" ˆ2 t u p ≠ " div(Ÿ 2 • G • (Òu p + ◊ p )) = f (2.1a) fl " 3 12 ˆ2 t ◊ p ≠ " 3 12 Div(C Á(◊ p )) + " Ÿ 2 • G • (Òu p + ◊ p ) = µ (2.1b) u p = ◊ p = 0 on ˆ (Clamped) (2.1c)
and for all symmetrical strain tensor Á

C Á = Q c c c c c a E x 1 ≠ ‹ xy ‹ yx ≠ E y ‹ xy 1 ≠ ‹ xy ‹ yx 0 ≠ E x ‹ yx 1 ≠ ‹ xy ‹ yx E y 1 ≠ ‹ xy ‹ yx 0 0 0 2 G xy R d d d d d b Q c a Á xx Á yy Á xy R d b , G = A G xz G yz B , Ÿ 2 = A Ÿ 2 x Ÿ 2 y B f and µ = A µ 1 µ 2 B
are the source terms of flexion and sheer.

Model 2.1 (Reissner-Mindlin) To condense the writings, we rewrite this system in the following matrix form

I M ¨ p + R p = F p in (2.2a) p = 0 on ˆ (2.2b) where p = A u p ◊ p B , M = A fl" 0 0 fl " 3 12 B , R = A B C ú C A B , F p = A f µ B (2.3)

Inria

The operators A, B and C are defined as

Y _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ [ A◊ = ≠ " 3 12 Div(CÁ(◊)) + "Ÿ 2 • G • ◊ (2.4a) Bu = ≠" div(Ÿ 2 • G • Òu) (2.4b) Cu = " Ÿ 2 • G • Òu (2.4c) C ú ◊ = ≠" Ÿ 2 • G • div ◊ (2.4d)

Modal decomposition

The R operator involving second-order derivatives in space is self-adjoint and can therefore be diagonalized.

Let D R = Ó p oe ! H 1 ( ) " 3 - -R p oe L 2 ( ) 3
Ô be the domain of the operator R.

As detailed in [Chabassier, 2012, p. 80], the spectral theorem ensures that there exists an orthonormal sequence (w n ) noeN ú of elements of D R which forms a Hilbertian basis of ! L 2 ( ) " 3 and a sequence

(⁄ n ) noeN ú of positive real numbers which tends to +OE such that 'n oe N ú , Rw n = ⁄ n Mw n (2.5)
The solution p and the source term F can then be decomposed in the basis

(w n ) noeN ú '(x, y, t) oe ◊ R + , Y _ _ _ ] _ _ _ [ p (x, y, t) = ÿ noeN ú p,n (t)w n (x, y) (2.6a) F p (x, y, t) = ÿ noeN ú F p,n (t)w n (x, y) (2.6b)
which allows us to rewrite the Reissner-Mindlin equations as a system of decoupled ordinary di erential equations (ODEs) verified by the modal amplitudes

'n oe N ú , ¨ p,n + ⁄ n p,n = F p,n (2.7)
This description already shows the numerical interest of modal decomposition by using a truncated system of N m modes, and the method of resolution will consist of determining the N m modal amplitudes of the solution at each simulation time step.

We refer to [Chabassier, 2012, p. 82] for detailed explanations of the damping phenomena that take place in the soundboard. Just note that we use a modal damping model intrinsic to each mode, for which each frequency is attenuated di erently according to an experimental model based on measurements by [Ege, 2009] and [Bucur, 2006]. It takes the form f ve (⁄) = a⁄ + b Ô ⁄ + c and allows us to rewrite the equations (2.7) with a viscoelastic damping factor

'n oe N ú , ¨ p,n + f ve (⁄ n ) ˙ p,n + ⁄ n p,n = F p,n (2.8)
More generally, a linear operator f ve (R) can be defined such that for any eigenmode

w n , f ve (R)w n = f ve (⁄ n )Mw n .
This gives a damped matrix formulation of the Reissner-Mindlin system:

Seek p : ◊ [0, T ] ≠ae R 3 so that I M ¨ p + f ve (R) ˙ p + R p = F p in (2.9a) p = 0 on ˆ (2.9b)
Model 2.2 (Abstract damped Reissner-Mindlin formulation)

Variational formulation and energy

Let us introduce the variational space

L D = ; p oe 1 H 1 ( ) 2 3 - - - -'k oe J D , ( p ) k = 0 on ˆ <
The set J D is used to specify which components are subject to Dirichlet conditions. For clamped edge conditions, J D = {1, 2, 3} and the space is simply

L D = ! H 1 0 ( ) " 3 .
For the simply supported conditions we'll use for simulations, J D = {1} which allows us to impose the Dirichlet condition only on the u p elevation motion, the other components being subject to Neumann conditions.

Seek p oe L D so that for all ú p oe L D ⁄ M ¨ p • ú p + ⁄ f ve (R) ˙ p • ú p + ⁄ R p • ú p = ⁄ F p • ú p (2.10)
In accordance with the previous paragraph, all modes (⁄ n , w n ) noeN ú satisfy at the variational level

'n oe N ú , ' ú p oe L D , ⁄ Rw n • ú p = ⁄ n ⁄ Mw n • ú p (2.11) Any regular enough solution p oe C 1 ([0, T ], L D ) fl C 2 ! [0, T ], L 2 ([0, L]) 3 " to the damped Reissner-Mindlin model 2.2 verifies Y _ _ _ ] _ _ _ [ dE p dt = ⁄ F p • ˙ p ≠ . . . ˙ p . . . 2 fve(R) (2.12a) E p (t) = 1 2 . . . ˙ p . . . 2 M + 1 2 Î p Î 2 R (2.12b)
Theorem 2.1 (Power balance of the damped Reissner-Mindlin system)

Proof. We use ú p = ˙ p oe L D in the variational formulation (2.10), and the result follows directly from this. The existence of a solution p oe C

1 ([0, T ], L D ) fl C 2 ! [0, T ], L 2 ( ) 3 "
follows from the application of the Hille-Yosida theorem to the Reissner-Mindlin system. The full proof is given in [Grob, 2006].

Bridge model with rotation mobilities

In [Chabassier, 2012], the choice was made to model only the vertical movement of the bridge (direction ‹ ) with a downbearing angle -. Simulation results show, however, that longitudinal modes are greatly under-represented in the sound produced by such a model.

We therefore propose a model that makes the bridge mobile around its point of contact C with the table. At this point, the moments of the forces exerted by the strings are transmitted to the shear directions ◊ 1 p and ◊ 2 p of the Reissner-Mindlin model. Indeed, the red bridge in figure 2, which has a height ¸typically between 38 and 64 millimeters, acts as a lever for these forces. For our modelling purposes, we assume that the bridge has no deformation, and therefore no energy, nor does it dissipate any. It is simply seen as a rigid body that acts as a coupling condition between the string and the soundboard. The reference frame (÷, • , ‹ ) is a Cartesian frame associated with the soundboard. Consequently, a displacement imposed on the bridge in the • direction is transmitted purely to the ◊ 1 shear angle, while a displacement on ÷ is transmitted purely to ◊ 2 . We therefore add a lateral angleto take into account any mismatch between the direction of the strings and the preferred directions of the wood defined by the stress tensor C of the model 2.1 and the presence of sti eners.

These two anglesand -, which orient the string relatively to the soundboard plane, are such that

Q c a ‹ • ÷R d b = Q c a cos(-) s i n ( -) ≠ sin(-) cos(-) cos(-) cos(-) sin(-) sin(-) ≠ cos(-) sin(-) R d b A ęu ęv B (3.1)

Expression of the coupling forces

Let us now denote F P , G P and H P the projections of the force exerted by the strings on the bridge in the three spatial directions (‹ , • , ÷).

They are taken into account in the right-hand-side of the Reissner-Mindlin equations (2.2) and distributed over the table by a distribution function ‰ :

M ¨ p + f ve (R) ˙ p + R p = F p = Q c a ≠F P (t)‰ (x, y) ≠¸G P (t)‰ (x, y) ≠¸H P (t)‰ (x, y) R d b (3.2)
To guarantee stable couplings, we impose an exact conservation of energy: the work transmitted by the strings must be equal to that injected into the soundboard. Thanks to the theorems A.1 and 2.1, it writes

T L • q(x = L, t) + ⁄ F p • ˙ p = 0 (3.3) with T L = Ë Aˆxq + Bq + ÒU(ˆxq) + Hˆx qÈ (x = L, t)
This is achieved by imposing continuity of the speeds at the bridge:

Y _ _ _ _ _ _ ] _ _ _ _ _ _ [ q(L, t) • ‹ c = ⁄ up ‰ (3.4a) q(L, t) • • c = ¸⁄ ◊1 p ‰ (3.4b) q(L, t) • ÷ c = ¸⁄ ◊2 p ‰ (3.4c)
where ‹ c , • c and ÷ c are the vectors used to select those components of the string unknown that have physical positional significance. For example, for the Timoshenko system, we have

q(L, t) • ‹ c = Q c a uc (L, t) 0 0 R d b • ‹ (3.5)
The physically dimensionless shear component Ï c , which is not involved in calculating the end-point velocity, has no projection in this base. For the geometrically exact model we have

q(L, t) • ‹ c = Q c a uc (L, t) vc (L, t) 0 R d b • ‹ (3.6)
The energy conservation equation (3.3) can then be rewritten as follows

T L • q(L, t) = ≠ ⁄ F p • ˙ p (3.7) = F P (t) ⁄ up ‰ + ¸GP (t) ⁄ ◊1 p ‰ + ¸HP (t) ⁄ ◊2 p ‰ (3.8) = [ q(L, t) • ‹ c ] F P (t) + [ q(L, t) • • c ] G P (t) + [ q(L, t) • ÷ c ] H P (t) (3.9)
By identification, the expression of the coupling forces are obtained

Y _ ] _ [ F P (t) = T L • ‹ c (3.10a) G P (t) = T L • • c (3.10b) H P (t) = T L • ÷ c (3.10c) Inria Seek q : [0, L] ◊ [0, T ] ≠ae R N , p : ◊ [0, T ] ≠ae R 3 so that Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [
String equation:

ˆ2 t Mq + ˆt(Rq ≠ ˆx(Hˆxq)) ≠ ˆx(Aˆxq + Bq + ÒU(ˆxq)) + Cq + t Bˆxq = F (x, t) (3.11a)
Boundary conditions for the string:

I 'k oe I D , q k (0, t) = 0 'k oe I N , ˆxq k (0, t) = 0 and Y _ _ _ _ _ _ ] _ _ _ _ _ _ [ q(L, t) • • c = ⁄ up ‰ q(L, t) • ‹ c = ¸⁄ ◊1 p ‰ q(L, t) • ÷ c = ¸⁄ ◊2 p ‰ (3.11b)
Coupling forces exerted on the soundboard :

Y _ _ ] _ _ [ F P (t) = T L • ‹ c G P (t) = T L • • c H P (t) = T L • ÷ c où T L = Ë Aˆxq + Bq + ÒU(ˆxq) + Hˆx qÈ (L, t) (3.11c)
Soundboard equation:

M ¨ p + f ve (R) ˙ p + R p = F p = Q c a ≠F P (t)‰ (x, y) ≠¸G P (t)‰ (x, y) ≠¸H P (t)‰ (x, y) R d b in (3.11d)
Boundary conditions for the soundboard:

'k oe J D , ( p ) k = 0 on ˆ (3.11e) Model 3.1 (String -Soundboard coupling)
This model now allows bridge mobility in the ‹ c and ÷ c directions, which was not the case in the work of [Chabassier, 2012]. It also allows coupling to a soundboard, which was not the case in [Tan, 2017]. The coupling is achieved by imposing a continuity of velocities between the end of the string attached to the bridge and the soundboard, which also ensures conservation of energy.

Variational formulation and energy

Let us introduce the variational spaces involved

Q D = ; q oe 1 H 1 ([0, L]) 2 N - - - -'k oe I D , q k (0) = 0 <
for the string

L D = ; p oe 1 H 1 ( ) 2 3 - - - -'k oe J D , ( p ) k = 0 on ˆ < for the soundboard RR n°9530 Seek q : [0, T ] ≠ae Q D , p : [0, T ] ≠ae L D and F P , G P , H P : [0, T ] ≠ae R so that for all (q ú , ú p ) oe Q D ◊ L D : Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [
String equation:

⁄ L 0 M q • q ú + ⁄ L 0 R q • q ú + ⁄ L 0 Hˆx q • ˆxq ú + ⁄ L 0 A C t B B A B A q ˆxq B • A q ú ˆxq ú B + ⁄ L 0 ÒU(ˆxq) • ˆxq ú = ⁄ L 0 F • q ú + F P (t) Ë q ú (L, t) • ‹ c È + G P (t) Ë q ú (L, t) • • c È + H P (t) Ë q ú (L, t) • ÷ c È (3.12a)
Coupling equations:

Y _ _ _ _ _ _ ] _ _ _ _ _ _ [ q(L, t) • ‹ c = ⁄ up ‰ q(L, t) • • c = ¸⁄ ◊1 p ‰ q(L, t) • ÷ c = ¸⁄ ◊2 p ‰ (3.12b) Soundboard equation: ⁄ M ¨ p • ú p + ⁄ f ve (R) ˙ p • ú p + ⁄ R p • ú p = ≠ F P (t) ⁄ u ú p ‰ ≠ G P (t)¸⁄ ◊ 1ú p ‰ ≠ H P (t)¸⁄ ◊ 2ú p ‰ (3.12c)
The forces F P , G P and H P are Lagrange multipliers associated with coupling constraints in the three directions of space.

Remark 3.1

Consider a string simply polarized in the vertical direction ęu . If -© 0[fi], there is no projection of motion in the ÷ direction and the third coupling equation in the ÷ c direction boils down to 0 = 0, leading to an ill-posed system. To overcome this problem, the multiplier H P and the associated coupling equation on ÷ c are eliminated for simulations in which -© 0 [fi].

Any regular enough solution

q oe C 1 ([0, T ], Q D ) fl C 2 1 [0, T ], L 2 ([0, L]) N 2 and p oe C 1 ([0, T ], L D ) fl C 2 ! [0, T ], L 2 ( ) 3 " to the model 3.1 verifies Y _ ] _ [ dE c,p dt = ⁄ L 0 F • q ≠ ⁄ L 0 R q • q ≠ ⁄ L 0 Hˆx q • ˆx q ≠ . . . ˙ p . . . 2 fve(R) (3.13a) E c,p (t) = E c (t) + E p (t) (3.13b)
where E c is the energy of the string defined by theorem A.1 and E p is the energy of the soundboard defined by theorem 2.1.

Theorem 3.1 (Power balance of the string -soundboard system)

Proof. The power balance of the string A.1 ensures that

dE c dt = T L • q(L, t) + ⁄ L 0 F • q ≠ ⁄ L 0 R q • q ≠ ⁄ L 0 Hˆx q • ˆx q (3.14)
and that of the soundboard 2.1 that

dE p dt = ⁄ F p • ˙ p ≠ . . . ˙ p . . . 2 fve(R) (3.15) = ≠F P (t) ⁄ up ‰ ≠ G P (t)¸⁄ ◊1 p ‰ ≠ H P (t)¸⁄ ◊2 p ‰ ≠ . . . ˙ p . . . 2 fve(R) (3.16)
According to the definitions of coupling forces (3.11d) we have

T L • q(L, t) = F P (t) Ë q(L, t) • ‹ c È + G P (t) Ë q(L, t) • • c È + H P (t) Ë q(L, t) • ÷ c È (3.17)
then the coupling conditions (3.12b) allow to write that

T L • q(L, t) = F P (t) ⁄ up ‰ + G P (t)¸⁄ ◊1 p ‰ + H P (t)¸⁄ ◊2 p ‰ (3.18)
The sum of the two power balances concludes the demonstration. By construction, all coupling terms vanish to guarantee a coupling that transfers energy exactly between the two subsystems.

Choice of distribution function ‰

The ‰ distribution function is used to distribute the force exerted by the string on the soundboard, and thus to convert a point force into a surface force. Several choices are possible, but all must verify that the integral of ‰ is 1 over R 2 , in order to achieve a weighted average.

In [Chabassier, 2012], ‰ is a Gaussian function defined by

'(x, y) oe , ‰ (x, y) = " fir 2 0 e ≠" ! r r 0 " 2 where r = . . . . . A x ≠ x 0 y ≠ y 0 B. . . . . 2 (3.19)
where " is a rate of decay, (x 0 , y 0 ) the center of the distribution, and r 0 the radius at which the influence becomes negligible. This choice is not, however, the most physically suitable for modelling a bridge, which is not very deformable, as the decay of the Gaussian is rather slow and doesn't allow the distribution zone to be focused precisely.

To overcome this shortcoming, we use another normalized distribution function, but one whose decay slope can be adjusted much more precisely

'(x, y) oe , Y _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ [ A x ỹB = A cos -‰ ≠ sin -‰ sin -‰ cos -‰ B A x y B (3.20a) " r,s (x) = 1 2r 3 1 1 + e ≠s(x+r) ≠ 1 1 + e ≠s(x≠r) 4 (3.20b) ‰ (x, y) = " rx,sx (x)" ry,sy (ỹ) (3.20c)
The parameters s x and s y are used to adjust the decay rates in the two directions to model the bridge's non-symmetry, and to orient it on the soundboard at an angle -‰ . In the bridge's width, the decay is very fast, and much slower in its length.

r x [m] r y [m] s x [-] s y [-] 0,02 0,05 1000 100 x 0 [m] y 0 [m] -‰ [°] 0 0 60 Figure 3: Distribution function ‰ .

Discretization of the coupled string -soundboard model

The discretization techniques used in this section are taken from [Derveaux, 2002] and [Chabassier, 2012]. They have been adapted to take the bridge rotation mobilities into account.

Numerical scheme for the soundboard

4.1.a Space discretization -Modal method

Let us consider h a conformal quadrangular mesh of and note

L h = Ó h oe L D - - -'K oe h , h|K oe Q r [X, Y ] 3 Ô where Q r [X, Y ]
is the space of polynomials generated by

) X i Y j | 1 6 i, j 6 r * .
In order to make the mass matrix diagonal, Gauss-Lobatto points are used to position the n h degrees of freedom of the finite element method, as well as the quadrature points to evaluate the integral terms. The basis functions

(" i ) ioe[ [1,n h ] ] are such that for any vertex S j , " i (S j ) = " ij .
The modal problem (2.11) can then be written in discrete form: For all n oe N ú , seek w h,n oe L h and ⁄ h,n oe R so that for all ú h oe L h :

⁄ Rw h,n • ú h = ⁄ h,n ⁄ Mw h,n • ú h (4.1)
Evaluating this variational formulation for all basis functions allows us to rewrite this problem in matrix form

R h W h,n = ⁄ h,n M h W h,n (4.2)
where W h,n is the vector of components of w h,n in the finite element basis. Since the R h and M h matrices are symmetrical and diagonal positive-definite respectively, the R h matrix is diagonalizable in a M h -orthogonal basis by the spectral theorem. If we denote P h the transition matrix from the finite element basis to this new modal basis (W h,n ) noeN ú we have

Y _ ] _ [ t P h R h P h = D h = diag (⁄ h,n ) t P h f ve (R h )P h = A h = diag (f ve (⁄ h,n )) t P h M h P h = I
and any h oe L h element can be expressed in the finite element basis EF h or in the modal basis mod h by the transformations

I EF h = P h mod h mod h = t P h M h EF h
The variational formulation of the damped Reissner-Mindlin model (2.10) can then be written in matrix form

M h ¨ EF h + f ve (R h ) ˙ EF h + R h EF h = F h (4.5)
and then in the modal basis

¨ mod h + A h ˙ mod h + D h mod h = t P h F h (4.6)
The modal resolution method then consists in truncating the sequence of modes to retain only the first N m , thus considerably reducing the calculation time compared with the number of degrees of freedom of the finite element basis. At the same time, very high-frequency finite element modes, which are polluted by discretization errors anyway, are eliminated [Durufle, 2006].

4.1.b Analytical resolution scheme in time

The time scheme we use consists in considering the second member constant from one time step to the next, which is a legitimate approximation since it is itself derived from a discretization of the strings and therefore no longer depends continuously on time.

Then let us consider that

F h (t) = F h (t n ) = F n h for all t oe [t n≠1/2 , t n+1/2 ].
At each and every iteration, seek

mod,n+1/2 h = mod h (t n+1/2 ) oe R n h so that I 't oe [t n≠1/2 , t n+1/2 ], ¨ mod h (t) + A h ˙ mod h (t) + D h mod h (t) = t P h F n h (4.7a) mod h (t n≠1/2 ) = mod,n≠1/2 h , ˙ mod h (t n≠1/2 ) = ˙ mod,n≠1/2 h (4.7b)
Numerical scheme 4.1 (Scheme for the soundboard)

It is in fact a system of decoupled ordinary di erential equations whose solution at time t n+1/2 is given analytically

(see appendix C) by Y ] [ mod,n+1/2 h = mod h (t n+1/2 ) = S 0 t mod h (t n≠1/2 ) + S 1 t ˙ mod h (t n≠1/2 ) + R t t P h F n h (4.8a) ˙ mod,n+1/2 h = ˙ mod h (t n+1/2 ) = Ṡ0 t mod h (t n≠1/2 ) + Ṡ1 t ˙ mod h (t n≠1/2 ) + Ṙ t t P h F n h (4.8b)
The solution to the scheme 4.1 verifies

Y _ _ _ ] _ _ _ [ E n+1/2 p ≠ E n≠1/2 p t = t P h F n h • mod,n+1/2 h ≠ mod,n≠1/2 h t ≠ 1 t ⁄ t n+1/2 t n≠1/2 . . . ˙ mod h . . . 2 A h (4.9a) E n+1/2 p = 1 2 . . . ˙ mod,n+1/2 h . . . 2 2 + 1 2 . . . mod,n+1/2 h . . . 2 D h (4.9b)
Theorem 4.1 (Discrete power balance of the soundboard modal system)

Proof. Take the scalar product of equation (4.7a) with ˙ mod h (t) and integrate for t oe [t n≠1/2 , t n+1/2 ]:

⁄ t n+1/2 t n≠1/2 ¨ mod h • ˙ mod h + ⁄ t n+1/2 t n≠1/2 D h mod h • ˙ mod h = ⁄ t n+1/2 t n≠1/2 t P h F n h • ˙ mod h ≠ ⁄ t n+1/2 t n≠1/2 A h ˙ mod h • ˙ mod h (4.10)
After integration, and since F h is constant on the considered interval, the equation reads:

1 2 3 . . . ˙ mod,n+1/2 h . . . 2 2 ≠ . . . ˙ mod,n≠1/2 h . . . 2 2 4 + 1 2 3 . . . ˙ mod,n+1/2 h . . . 2 D h ≠ . . . ˙ mod,n≠1/2 h . . . 2 2 4 = t P h F n h • 1 mod,n+1/2 h ≠ mod,n≠1/2 h 2 ≠ ⁄ t n+1/2 t n≠1/2 . . . ˙ mod h . . . 2 A h (4.11)
The result of the theorem is obtained by simply dividing by t.

Inria

Numerical scheme for the coupled model

We can now deal with the coupling between the string and the soundboard. In order to remain generic, the numerical scheme used for the string is denoted C. It is assumed to be centred at time t n , for example 2-IEQ, 2-SAV [START_REF] Castera | Numerical analysis of quadratized schemes. Application to the simulation of the nonlinear piano string[END_REF] or the discrete gradient schemes developed in [Chabassier, 2012]. The use of the P-SAV scheme proposed in [START_REF] Castera | Numerical analysis of quadratized schemes. Application to the simulation of the nonlinear piano string[END_REF] requires slight adaptation of the schemes to ensure good consistency in t n+1/2 by discretizing all variables, including those of the soundboard, onto a single {t n } grid. The 2-SAV scheme that will be used for the following simulations is presented in the appendix A.

Let's start by defining discrete analogues of the positioning vectors ‹ c , • c and ÷ c , which we call ‹ c,h ,

• c,h and ÷ c,h such that

'q ú h oe Q h , Y _ _ ] _ _ [ q ú h (L, t) • ‹ c = Q ú h • ‹ c,h q ú h (L, t) • • c = Q ú h • • c,h q ú h (L, t) • ÷ c = Q ú h • ÷ c,h (4.12)
and where Q ú h is the vector of components of q ú h in the finite element basis of the string.

J h is the vector that distributes the coupling on the soundboard in a discrete manner

'i oe [[1, n h ]], (J h ) i = j h ‰ (x, y)" i (x, y)dxdy (4.13) so that EF h • Q c a J h 0 0 R d b = j h u p,h ‰ EF h • Q c a 0 J h 0 R d b = j h ◊ 1 p,h ‰ EF h • Q c a 0 0 J h R d b = j h ◊ 2 p,h ‰ (4.14)
With these definitions, the velocity continuity conditions (3.4) can be written at the discrete level

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ Qh • ‹ c,h = ˙ EF h • Q c a J h 0 0 R d b = ˙ mod h • t P h Q c a J h 0 0 R d b = ˙ mod h • t P U h J h (4.15a) Qh • • c,h = ¸˙ EF h • Q c a 0 J h 0 R d b = ¸˙ mod h • t P h Q c a 0 J h 0 R d b = ˙ mod h • ¸tP 1 h J h (4.15b) Qh • ÷ c,h = ¸˙ EF h • Q c a 0 0 J h R d b = ¸˙ mod h • t P h Q c a 0 0 J h R d b = ˙ mod h • ¸tP 2 h J h (4.15c)
where P U h , P 1 h and P 2 h are the transition submatrices associated with the di erent components U p,h , 1 p,h and 2 p,h of the solution vector h

P h = Q c c c c c c c c c a P U h P 1 h P 2 h R d d d d d d d d d b t P h = Q c a t P U h t P 1 h t P 2 h R d b
At each and every iteration, seek

Q n+1 h oe R n q h , mod,n+1/2 h oe R n h and (F P , G P , H P ) oe R 3 so that Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [ String scheme: C(Q n+1 h , Q n h , Q n≠1 h ) = F n h + F P ‹ c,h + G P • c,h + H P ÷ c,h (4.16a)
Coupling scheme:

Y _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ [ Q n+1 h ≠ Q n≠1 h 2 t • ‹ c,h = mod,n+1/2 h ≠ mod,n≠1/2 h t • t P U h J h Q n+1 h ≠ Q n≠1 h 2 t • • c,h = mod,n+1/2 h ≠ mod,n≠1/2 h t • ¸tP 1 h J h Q n+1 h ≠ Q n≠1 h 2 t • ÷ c,h = mod,n+1/2 h ≠ mod,n≠1/2 h t • ¸tP 2 h J h (4.16b)
Soundboard scheme:

Y _ _ _ ] _ _ _ [ 't oe [t n≠1/2 , t n+1/2 ], ¨ mod h (t) + A h ˙ mod h (t) + D h mod h (t) = t P h F n h = ≠F P t P U h J h ≠ G P ¸tP 1 h J h ≠ H P ¸tP 2 h J h mod h (t n≠1/2 ) = mod,n≠1/2 h , ˙ mod h (t n≠1/2 ) = ˙ mod,n≠1/2 h (4.16c)

Numerical scheme 4.2 (Coupling string -soundboard)

The solution to the scheme 4.2 verifies

Y _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ [ E n+1/2 c,p ≠ E n≠1/2 c,p t = F n h • Q n+1 h ≠ Q n≠1 h 2 t ≠ . . . . . Q n+1 h ≠ Q n≠1 h 2 t . . . . . 2 R h ≠ 1 t ⁄ t n+1/2 t n≠1/2 . . . ˙ mod h . . . 2 A h (4.17a) E n+1/2 c,p = E n+1/2 c + E n+1/2 p (4.17b)
where E c is the discrete scheme-dependent energy of the string. See for example the theorem A.3 for the 2-SAV scheme. E p is the soundboard energy defined by the theorem 4.1.

Theorem 4.2 (Discrete power balance of the scheme 4.2)

Proof. The proof is similar to that of the continuous theorem 4.2: all coupling terms simplify to guarantee exact coupling in terms of energy exchange.

Inria

Solving algorithm

The 4.2 scheme discretizes the coupling equations implicitly.

That being said, it is possible to make the resolution method explicit by using the equation (4.8a):

Y _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ [ Q n+1 h ≠ Q n≠1 h 2 • ‹ c,h = 1 S 0 t mod,n≠1/2 h + S 1 t ˙ mod,n≠1/2 h + R t t P h F n h ≠ mod,n≠1/2 h 2 • t P U h J h Q n+1 h ≠ Q n≠1 h 2 • • c,h = 1 S 0 t mod,n≠1/2 h + S 1 t ˙ mod,n≠1/2 h + R t t P h F n h ≠ mod,n≠1/2 h 2 • ¸tP 1 h J h Q n+1 h ≠ Q n≠1 h 2 • ÷ c,h = 1 S 0 t mod,n≠1/2 h + S 1 t ˙ mod,n≠1/2 h + R t t P h F n h ≠ mod,n≠1/2 h 2 • ¸tP 2 h J h
Let us introduce

F n≠1/2 = S 0 t mod,n≠1/2 h + S 1 t ˙ mod,n≠1/2 h (4.18)
which allows the coupling conditions to be rewritten in explicit form:

Y _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ [ Q n+1 h ≠ Q n≠1 h 2 • ‹ c,h = A n≠1/2 1 ≠ F P B 1,1 ≠ G P B 1,2 ≠ H P B 1,3 (4.19a) Q n+1 h ≠ Q n≠1 h 2 • • c,h = A n≠1/2 2 ≠ F P B 2,1 ≠ G P B 2,2 ≠ H P B 2,3 (4.19b) Q n+1 h ≠ Q n≠1 h 2 • ÷ c,h = A n≠1/2 3 ≠ F P B 3,1 ≠ G P B 3,2 ≠ H P B 3,3 (4.19c)
where A n≠1/2 is the vector defined by

A n≠1/2 = Q c c c c c c c a 1 F n≠1/2 ≠ mod,n≠1/2 h 2 • t P U h J h 1 F n≠1/2 ≠ mod,n≠1/2 h 2 • ¸tP 1 h J h 1 F n≠1/2 ≠ mod,n≠1/2 h 2 • ¸tP 2 h J h R d d d d d d d b (4.20)
and B the constant matrix

B = R t Q c c c c c a t P U h J h • t P U h J h ¸tP 1 h J h • t P U h J h ¸tP 2 h J h • t P U h J h t P U h J h • ¸tP 1 h J h ¸tP 1 h J h • ¸tP 1 h J h ¸tP 2 h J h • ¸tP 1 h J h t P U h J h • ¸tP 2 h J h ¸tP 1 h J h • ¸tP 2 h J h ¸tP 2 h J h • ¸tP 2 h J h R d d d d d b (4.21)
The method for solving the 4.2 scheme, which treats the string and soundboard decoupled, is then summarized by the following algorithm:

Y _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ [
Initialize all the variables and pre-compute S 0 t , S 1 t , R t and B For each and every iteration:

-Update F n≠1/2 and A n≠1/2 -Solve the string's scheme with the coupling conditions:

Y _ _ _ _ _ _ _ _ _ _ _ ] _ _ _ _ _ _ _ _ _ _ _ [ C(Q n+1 h , Q n h , Q n≠1 h ) = S n h + F P ‹ c,h + G P • c,h + H P ÷ c,h Q n+1 h ≠ Q n≠1 h 2 • ‹ c,h = A n≠1/2 1 ≠ F P B 1,1 ≠ G P B 1,2 ≠ H P B 1,3 Q n+1 h ≠ Q n≠1 h 2 • • c,h = A n≠1/2 2 ≠ F P B 2,1 ≠ G P B 2,2 ≠ H P B 2,3 Q n+1 h ≠ Q n≠1 h 2 • ÷ c,h = A n≠1/2 3 ≠ F P B 3,1 ≠ G P B 3,2 ≠ H P B 3,3 (4.22a)
-Update the soundboard unknowns:

mod,n+1/2 h = F n≠1/2 ≠ F P R t t P U h J h ≠ G P R t ¸tP 1 h J h ≠ H P R t ¸tP 2 h J h (4.22b)

Numerical results

In this section, we present comparisons between the refined model with three degrees of freedom in the ‹ , • and ÷ directions, hereafter referred to as "3dof", and the [Chabassier, 2012] model where the bridge has only one degree of freedom in the ‹ direction, hereafter referred to as "1dof". Note that theangle has no influence on the 1dof model, for which it is always set to 0, and the

H P multiplier is eliminated.
After a paragraph on calculating the soundboard modes, we check the energy conservation and convergence properties of the new scheme.

The parameters used for the following simulations are summarized below. In place of the hammer, the geometrically exact and sti ness-free string F 3 is excited by a C OE source in space and time as defined by appendix B.

String's scheme

Scheme N x ◊ Order Auxiliary constant Stabilization Damping 2-SAV 11 1/4 4 10 4 -ú No Source term A [-] x 0 [m] ‡ x [m] t 0 [ms] ‡ t [ms]
10000 0,115 0,01 2,5 1,5

Inria

The soundboard is a 9-millimeter-thick spruce plate whose wood fibers are oriented at an angle of -= 54.4°. 16 spruce sti eners, also oriented at -= ≠40°and 25 millimeters thick, are glued to it.

They are shown in green of the adjacent drawing.

The pink bridge of 40 mm high, is made of beech. These thickness variations are taken into account by locally changing the thickness of the Reissner-Mindlin plate.

Wood fl

E xx E yy ‹ xy ‹ yx G xy G xz G yz [kg.m ≠3 ] [GPa] [GPa] [-] [-] [GPa] [GPa] [GPa]
Spruce 380 11,0 0,65 0,26 0,015 0,66 1,2 0,042 Beech 750 13,7 2,24 0,3 0 1,61 1,06 0,46 The parameters used to model the soundboard damping in simulations where it is taken into account are such that (see paragraph 2.2) the modes of the soundboard are computed by the Arpack1 software from a 2D mesh and with 4-th order finite elements.

f ve (f ) = 2 ◊ 10 ≠5 f 2 + 7 ◊ 10 ≠2 f (4.23)

4.4.a Computation of the modes

Although the mesh in figure 4 is composed of triangles and quadrangles, the Montjoie2 finite element code automatically cuts triangles into quadrangles to obtain a purely quadrangular mesh and allow mass lumping. This increases the solver's e ciency, as the mass matrix M h doesn't have to be factorized. In figure 4, the sti eners are shown in green, and the bridge in pink.

Some modal deformations from simulations are shown in figure 5. They are consistent with the Chladni shapes obtained experimentally in [Conklin, 1996] as well as those measured and simulated in step 3 of [START_REF] Corradi | Modal analysis of a grand piano soundboard at successive manufacturing stages[END_REF].

The influence of the sti eners and bridge can be seen in figure 5 in the form of hatched patterns clearly visible on certain angular modes (modes 2, 3 and 7).

Finally, it should be noted that the computation of some 2,500 modes is required to reach the 10 kHz frequency.

4.4.b Numerical preservation of the power balance

In the absence of damping and for -= 45°and -= 45°, we define the energy residual Á in relation to the theorem 4.2 as

Á = E n+1/2 c,p ≠ E n≠1/2 c,p E max ≠ F n h • Q n+1 h ≠ Q n≠1 h 2E max (4.24)
Figure 6: Energy residual for 1dof (left) and 3dof (right) transmission.

It is clear from figure 6 that the Á energy residuals are multiples of the machine error. The discrete power balance is therefore verified. Figure 7 shows time convergence curves of the 4.2 scheme for the 1dof transmission model and the 3dof refinement. Errors are calculated in ¸OE norm by dividing the time step consecutively by 2 and at specific positions on the string (x = L/4 and x = L) and on the table (P 2 position).

4.4.c Time convergence of the scheme

'(X, Y ) oe R n ◊ R n , E ¸OE (X, Y ) = ÎX ≠ Y Î OE ÎY Î OE = max ioe[ [1,n] ] - -x i ≠ y i - - max ioe[ [1,n] ] |y i | (4.25) Inria
Apart from the first point with a large time step, which is heavily polluted due to the 2-SAV string scheme as described in [START_REF] Castera | Numerical analysis of quadratized schemes. Application to the simulation of the nonlinear piano string[END_REF], the convergence is quadratic. An increase in error, also described in [START_REF] Castera | Numerical analysis of quadratized schemes. Application to the simulation of the nonlinear piano string[END_REF] and typical of the 2-SAV method, is visible for small time steps.

Transmission of longitudinal waves and precursors

Let us remember that transverse and longitudinal waves have very di erent propagation velocities in the strings. Longitudinal waves usually travel about 15 times faster than transverse waves, giving rise to the longitudinal precursor phenomenon observed in [START_REF] Podlesak | Dispersion of waves in piano strings[END_REF].

The greater the di erence between the arrival times of the two waves at the bridge, the more visible the precursor.

In this section, we keep on considering a non-linear string with no sti ness, excited by a regular source that produces a fortissimo nuance. Sti ness is not taken into account in order to observe only longitudinal phenomena without polluting them with dispersive phenomena. We show results on a low D 1 string, a medium F 3 string about half as long, and an even shorter high C 5 string. We add damping in the strings and in the soundboard, and the string mesh uses N x = 70 nodes.

D 1 note

Inria

Figure 8 shows the di erent positions taken by the x = L end of the string. On the left, for a 1dof transmission, the string is constrained vertically along the ‹ axis. With 3dof, the movement is no longer strictly vertical, and a longitudinal component appears. In the figure on the right, with an angle -= 45°, we can make the same observations, but in a frame rotated by 45°, consistent with the string's orientation.

Figure 9 shows the acceleration of the x = L end of the string attached to the bridge in its transverse direction on the left and its longitudinal direction on the right for an angle -= 0. With a 1dof transmission, everything happens as if the string were subject to a Dirichlet condition in the longitudinal direction. As it can be seen in the figure on the right, the x = L end of the string has no mobility in the longitudinal direction, and only the transverse component visible in the figure on the left starting at 12 ms can be transmitted. With 3dof transmission, the string is free to move in its longitudinal direction, as it can be seen on the right-hand figure 9. Since longitudinal waves propagate much faster than transverse waves, the first longitudinal vibrations appear at 2 ms, i.e. much faster than the first transverse vibrations that were visible with 1dof. Consequently, the vibrations in the interval [2 ms, 12 ms] are entirely due to longitudinal phenomena, whether caused by the longitudinal waves themselves in the figure on the right, or by their non-linear coupling with the transverse waves in the figure on the left. The vibratory signal between 2 ms and 12 ms is a longitudinal precursor.

Precursor Figure 10: Soundboard acceleration at P 1 (left). Listening signal (right).

The same phenomena are visible on the soundboard and illustrated in figure 10. The figure on the left shows the acceleration of point P 1 on the soundboard, and that on the right a listening signal reconstructed by combining several points on the soundboard and recorded at a point distant from the piano (see appendix D). With a 1dof transmission, the first waves are visible from 12 ms at point P 1 and delayed to 17 ms at the listening point due to its distance. With 3dof transmission, waves are visible from 2 ms at point P 1 , and from 7 ms at the listening point. The additional vibrations constitute the longitudinal precursor that is potentially audible at the sound attack, given its high amplitude.

The power spectrum of the listening signal is shown in figure 11. Fairly significant amplitude di erences can be observed along the entire frequency range, but specifically marked on the longitudinal partials, which are influenced by the presence of the precursor. The observations made on the note D 1 remain valid, but the precursor is less visible and shorter (less than 2 ms) due to the reduced string length. The power spectrum of the listening signal in figure 15 shows fairly significant di erences, most pronounced on the longitudinal partials.

F 3 note

C 5 note

For this note, due to its short length (L = 32.6 cm), we use a less refined string mesh with N x = 20 nodes. Here again, the precursor exists, but has no oscillations of its own and lasts less than a millisecond, which is coherent with the short length of the string. The power spectrum of the listening signal is again quite strongly modified.

Figure 19: Power spectrum of the listening signal.

Inria

Summary

The model developed allows us to account for complex bridge movements in the soundboard elevation and shear directions, and also shows the presence of longitudinal precursors throughout the piano range. Their amplitude and duration are consistent with string length: the longer the strings, the longer the precursors and the higher their amplitude. Their importance decreases rapidly with string length.

Beyond the presence of precursors, there are also very significant di erences between signals with and without transmission. Spectral contents are significantly di erent, mainly in longitudinal resonances.

Influence of the string orientation relatively to the soundboard

Due to the presence of the oriented sti eners, the bridge and the anisotropic properties of wood, it is possible that the orientation of the strings relatively to the soundboard has an e ect on the transmission of waves at the bridge. This is studied in this section.

To do this, we study the admittance of the bridge, i.e. the possible mobility in the three free directions (‹ , ÷, • ) as a function of frequencies and as a function of the lateral angledefined by figure 2.

Bridge admittance -Frequency mobility

Let us go back to the equation that describes the soundboard's vibrations

M ¨ p + f ve (R) ˙ + R p = F (x, y, t) (6.1)
and apply a Fourier transform to study it in the frequency domain

1 ≠Ê 2 M + iÊf ve (R) + R 2 ˆ p = F (x, y, Ê) (6.
2)

The bridge admittance matrix [START_REF] Chaigne | Acoustique des instruments de musique[END_REF], [START_REF] Cuenca | Three-dimensional interaction between strings, bridge and soundboard in modern piano's treble range[END_REF] defined by

Ŷ(Ê) F (x, y, Ê) = iÊ ˆ p (x, y, Ê) (6.3) writes Ŷ(Ê) = iÊ 1 R + iÊf ve (R) ≠ Ê 2 M 2 ≠1 (6.4)
The same reasoning on the discrete modal problem gives an admittance matrix

Ŷh (Ê) = iÊ 1 D h + iÊA h ≠ Ê 2 I h 2 ≠1 (6.5)
For each pulse Ê, this matrix characterizes the soundboard's mobility as a function of the input force applied to it. Of course, this input force depends on the relative orientation of the string and the soundboard.

Power transmitted to the soundboard

The instantaneous power transmitted through the bridge is given by

P(t) = ⁄ F (x, y, t) • ˙ p (x, y, t)dxdy (6.6)
and Parseval's formula ensures equivalence of norms between the time domain and the Fourier domain:

⁄ R + P(t)dt = ⁄ R P(t)dt = 1 2fi ⁄ R ⁄ F ú (x, y, Ê) • iÊ ˆ p (x, y, Ê)dxdydÊ = 1 2fi ⁄ R Î(Ê)dÊ (6.7)
where ú denotes the conjugate complex. Note that Î(Ê) = Îú (≠Ê) since P(t) is real.

From this relation, we can define a power transmitted frequency by frequency as

P(Ê) = 1 2fi 1 Î(Ê) + Îú (≠Ê) 2 = 1 fi Re 1 Î(Ê) 2 = 1 fi ⁄ Re 1 F ú (x, y, Ê) • iÊ ˆ p (x, y, Ê) 2 dxdy (6.8)
where Re denotes the real part. With the exception of one multiplicative constant, the formula used in [START_REF] Suzuki | Vibration and sound radiation of a piano soundboard[END_REF] is the same.

Inria

The right-hand-side vector of equation (6.2), which is the force applied to the bridge, is given by

F (x, y, Ê) = ≠‰ (x, y) Q c a F P ¸Ĝ P ¸Ĥ P R d b = ≠‰ (x, y) Q c a t ‹ c ¸t• c ¸t÷ c R d b TL (w) = ≠R --(x, y) TL (Ê) (6.9)
where

R --= ≠‰ (x, y) Q c a t ‹ c ¸t• c ¸t÷ c R d b (6.10)
is the matrix used to characterize the orientation of the string in relation to the soundboard. The transmitted power is then expressed as a function of admittance

P--(Ê) = 1 fi ⁄ Re 1 ŶR --TL • R --TL 2 dxdy (6.11)
At the discrete level, this writes

P--,h (Ê) = 1 fi Re 1 Ŷh R --,h TLh • R --,h TLh 2 (6.12) with R --,h = t P U h J h t ‹ c,h + ¸tP 1 h J h t • c,h + ¸tP 2 h J h t ÷ c,h (6.13)
and where the vectors ‹ c,h , • c,h and ÷ c,h depend on the orientations (-, -) via the relation (3.1). This relation (6.12) allows us to observe, frequency by frequency, the influence of the string's lateral orientation on the power transmitted to the soundboard.

Transverse impulse response

Let us choose TL = t (1, 0) which corresponds to an impulse in the transverse direction of the string.

We choose a downbearing angle -= 4°and the impulse is applied on the soundboard at the F 3 string attachment point. By varying the lateral angleof the strings (see figure 2), we can draw a power transmission map P--,h (Ê) shown in figure 20. The quantity plotted on figure 20 is the power multiplied by the pulsation Ê P--,h (Ê) for better visualization. We note that the lateral anglehas very little influence on the power transmission, as shown by the invariance of the color on the vertical axis.

Longitudinale impulse response

Let's now choose TL = t (0, 1) which corresponds to an impulse in the longitudinal direction of the string. We choose a downbearing angle -= 4°and the impulse is applied on the soundboard at the F 3 string attachment point. By varying the lateral angle of thestrings (see figure 2), we can draw power transmission maps P--,h (Ê) represented on figures 21 and 22. On these maps, dark areas correspond to high power transmissions, while light areas correspond to low power transmissions. Let us take an example of how to read figure 21. Let us imagine that the frequency f ƒ 222 Hz (mode 14) framed in orange on the figure is annoying to the ear, too present in the sound, and that a piano maker wishes to attenuate it. In this case, a relevant choice of string orientation would be oe [0, 40 ¶ ], since, as shown by the very light color of the line for this choice of angle, very little power is transmitted. However, this possibility is limited by the layout constraints of the piano's various components, which already take a lot of space. Therefore in practice it is not possible to adjust the lateral angle of each string independently.

Inria

If we consider a wider range of frequencies in figure 22, we can see very clearly that darker areas of high power transmission and lighter areas of lower transmitted power form clusters (marked in dotted orange in figure 22). In the mid-frequency range for example, from 1 to 4 kHz, the -> 0 zone is much darker and transmits much more power than the much lighter -< 0 zone. In the high-frequency range, from 7 to 10 kHz, an imposing dark spot of high power transmission is visible at anglesoe [≠60 ¶ , 40 ¶ ].

The transmission map in figure 22 does not reveal any optimal orientations that would uniformly limit or maximize the transmitted power over the entire audible frequency range [20 Hz -20 kHz]. However, there is a clear di erence in behaviour between low and high frequencies. These di erences could be explained by the orientations of the wood fibres of the soundboard (-= 54.4°), sti eners (-= ≠40°) and bridge (-= 23°), which influence shear strength and could have an influence on power transmission. In fact, it is known that wood strength is higher for torsional forces applied orthogonally to the fibres and, conversely, lower for forces applied parallel to the fibres. To explain the areas of high and low power transmission in figure 22, it would be necessary to be able to assign an importance to the resistance imposed by each element, as a function of frequency. It's possible, for example, that the wood fibres of the table are small enough to be unnoticeable at low frequencies, but have a significant importance at high frequencies. And vice versa, the sti eners and bridge, which are quite massive, may have a stronger influence on low frequencies than on high frequencies, or else impose a confinement of vibrations as shown in [Ege, 2009].

We do not go any further into this part, which would require a very detailed study on its own, but we do think it could lead to interesting studies in terms of instrument making.

Conclusion and Prospects

In this report, we have developed a bridge model in the form of a coupling condition that allows longitudinal string vibrations to be transmitted to the soundboard. The numerical method is based on a modal decomposition of the soundboard and an analytical time resolution on a truncated mode sequence. The coupling scheme preserves a discrete energy balance that ensures its stability, and converges quadratically in time.

For realistic string angular configurations (-small), longitudinal waves are well transmitted to the soundboard and precursors are observed. Finally, the last section studied the soundboard's impulse responses. For an impulse in the transverse direction of the string, we found no significant influence of the lateralangle on the power transmitted via the bridge. On the other hand, string orientations clearly appear to maximize or minimize the transmission of longitudinal waves over certain frequency ranges. Further studies may enable us to optimize instrument making at this level of the instrument.

Although it allows us to account for the expected physical phenomena, in particular precursors that seem realistic, this model still lacks experimental validation. Experiments would make it possible to validate the use of these coupling conditions as a reduced model of the bridge, or, on the contrary, to envisage more complex modelling, for example by considering the bridge in its entirety as a beam and coupling all its degrees of freedom with those of the soundboard. This problem, typically the coupling of a Timoshenko beam to a Reissner-Mindlin plate, in the time domain and with good energy properties, remains to our knowledge an open problem today.

As mentioned previously in [Chabassier, 2012], our model still does not take into account the downbearing exerted on the soundboard in the resting position. Piano builders are well aware, as we mentioned in the introduction to this report, that adjusting the pressure of the strings on the soundboard at rest is one of the most fundamental criteria for sound quality. A pre-stressed shell model, as in [START_REF] Mamou-Mani | Numerical simulation of a piano soundboard under downbearing[END_REF], would perhaps make it possible to better model the transmission of (longitudinal) waves at the bridge. Once again, we lack recent and precise experiments at this level of the instrument to validate or invalidate our models.

Finally, on almost all pianos, the bridges are covered with a layer of graphite that allows the string to slide longitudinally, particularly when tuning, to avoid breaking the bridges. It is therefore possible that the longitudinal waves at the end of the strings are not transmitted entirely to the bridge, but that some of them are also transmitted to the lengths of dead string located after the bridge, which would be interesting to model.

At each and every iteration, seek Q n+1 h oe R n q h and z n+1/2 h oe R so that

I M h " 2 Q n h + R h "µQ n h + K h {Q h } n ◊ + µz n h G(Q n h ) = F n h (A.15a) "z n h = G(Q n h ) • "µQ n h (A.15b) and Q 0 h = Q 1 h = 0 and z 1/2 h = Ô c.

Numerical scheme A.1 (2-SAV scheme)

This scheme is linaerly implicit and verifies the following discrete power balance.

The solution to the scheme A.1 verifies A is the chosen amplitude. The parameters used to plot the figure B23 are presented in the table B3.

Y _ ] _ [ "E n c = F n h • "µQ n h ≠ Î"µQ n h Î 2 R h (A.16a) E n+1/2 c =
A [-] x 0 [m] ‡ x [m] t 0 [s] ‡ t [

Inria

C Exact solution of the soundboard's time scheme

Solving the numerical soundboard scheme is equivalent to solving the following diagonal system of ordinary di erential equations with constant second member F for each time step:

U (t) = S 0 t≠t 0 U 0 + S 1 t≠t 0 U 1 + R t≠t 0 F (C.2a) U (t) = Ṡ0 t≠t 0 U 0 + Ṡ1 t≠t 0 U 1 + Ṙt≠t 0 F (C.2b)
where the diagonal matrices S 0 , S 1 and R are defined by

If |-i | < ⁄ i If |-i | > ⁄ i ! S 0
|-2 i ≠ ⁄ 2 i |
The proof can be found in [Derveaux, 2002] and in [Chabassier, 2012].

D Generation of the soundboard's listening signals

During the simulations, seismograms of the soundboard are calculated at a list of points (P i ) such that P 1 = (0.17, 1.49) P 2 = (0.30, 0.30) P 3 = (0.34, 1.00) P 4 = (0.59, 1.26) P 5 = (0.60, 0.60) P 6 = (0.60, 1.70) P 7 = (1.00, 0.17)

The listening point is located in a realistic position for a listener, at (2.00, 2.00, 0.60), i.e. in the hollow of the grand piano, and 60 centimetres above the soundboard.

To generate the listening signal, we simply sum the soundboard acceleration signals at the points P i , weighted by the distance d i separating them from the listening point, and with a corresponding time o set. Inria

S(t) =

Figure 1 :

 1 Figure 1: A piano bridge.

Figure 2 :

 2 Figure 2: Diagram of string -bridge -soundboard coupling

Figure 4 :

 4 Figure 4: 2D-mesh of the soundboard.

Figure 5 :

 5 Figure 5: Modal soundboard deformations. Membrane component (left) and shear angles (center, right).

Figure 7 :

 7 Figure 7: Time convergence for 1dof (left) and 3dof (right) transmission.

Figure 8 :

 8 Figure 8: Motion of the end of the string with -= 0 (left) and -= 45°(right).

Figure 9 :

 9 Figure 9: Acceleration of the end of the string for -= 0. Transverse (left), longitudinal (right).

Figure 11 :

 11 Figure 11: Power spectrum of the listening signal.

Figure 12 :

 12 Figure 12: Motion of the end of the string with -= 0 (left) and -= 45°(right).

Figure 13 :

 13 Figure 13: Acceleration of the end of the string for -= 0. Transverse (left), longitudinal (right).

Figure 14 :

 14 Figure 14: Soundboard acceleration at P 1 (left). Listening signal (right).

Figure 15 :

 15 Figure 15: Power spectrum of the listening signal.

Figure 16 :

 16 Figure 16: Motion of the end of the string with -= 0 (left) and -= 45°(right).

Figure 17 :

 17 Figure 17: Acceleration of the end of the string for -= 0. Transverse (left), longitudinal (right).

Figure 18 :

 18 Figure 18: Soundboard acceleration at P 1 (left). Listening signal (right).

Figure 20 :

 20 Figure 20: Transmitted power Ê P--,h (Ê) with respect to -(vertical axis) and frequency (horizontal axis).

Figure 21 :

 21 Figure 21: Transmitted power P--,h (Ê) function ofand of frequency [0 -500 Hz].

Figure 22 :

 22 Figure 22: Transmitted power P--,h (Ê) function ofand of frequency [0 -10 kHz].

Figure B23 :

 B23 Figure B23: View of the source S for a fixed time (left) and a fixed position (right).

  diag(2i ) and D = diag(⁄ 2 i )The solution and its time derivative write Y ] [

  m/s is the speed of sound in the air surrounding the piano.

Table 1 :

 1 Discretization parameters for the string.

Table 2 :

 2 Mechanical parameters for the soundboard.

  Theorem A.3 (Discrete power balance of the 2-SAV scheme A.1)B C OE source term in space and in timeThe following function S is C OE in space and in time if(x, t) oe [x 0 ≠ ‡ x , x 0 + ‡ x ] ◊ [t 0 ≠ ‡ t , t 0 + ‡ t ] 0 elseThe source is applied on the space interval[x 0 ≠ ‡ x , x 0 + ‡ x ] centered in x 0 during a time 2 ‡ t on the interval [t 0 ≠ ‡ t , t 0 + ‡ t ].

	Y _ _ _ _ ]	A e 1≠	1≠	1 ! x≠x 0 ‡x	" 2	e 1≠	1 ‡ t ! t≠t 0 1≠	" 2
	where Ê M S(x, t) = _ _ _ _ [						1 2	. . ."Q n+1/2 h	. . . 2 Â M h	+	1 2	. . .µQ n+1/2 h	. . . 2 K h	+	1 2	1	z	n+1/2 h	2 2	(A.16b)

h = M h + t 2 1 ◊ ≠ 1 4 2

K h is the modified mass matrix.

Table B3 :

 B3 Parameters of the source.

					s]
	1	0,35	0,08	0,2	0,15

RR n°9530

https://en.wikipedia.org/wiki/ARPACK

https://www.math.u-bordeaux.fr/~durufle/montjoie/ RR n°9530

Inria

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

A Modelling and numerical simulation of the strings

This appendix is a condensed extract of [START_REF] Castera | Numerical analysis of quadratized schemes. Application to the simulation of the nonlinear piano string[END_REF], in which more detailed explanations about numerical methods for the strings can be found.

A.1 String models

There are many models that can be used to describe strings. Their complexity varies according to the accuracy required and the di erent physical phenomena taken into account (sti ness, nonlinearity, double polarisation). Most of the existing one-dimensional models, with the exception of the Euler-Bernoulli model, can be written in the same unified matrix form. For the list of models below, we recall the meaning of the physical coe cients. rho is the density, S the cross-section, T 0 the tension at rest, E the Young's modulus and G the shear modulus, I the quadratic moment of section, and Ÿ the Timoshenko coe cient.

R u , R v , R Ï designate dry damping, and ÷ u , ÷ v , ÷ Ï viscous damping.

Vibrating string model (linear non-sti )

Geometrically exact model (nonlinear non-sti )

Geometrically exact sti model (nonlinear sti )

A.2 Variational formulation and energy

Let us introduce the variational space Q D of the functions that cancels at x = 0:

The x = L end of the string is not subject to a Dirichlet condition, and we denote T L the tension exerted there by the external coupling. Seek q oe Q D so that for all q ú oe Q D :

Theorem A.1 (Power balance of the string model)

Proof. Just apply the variational form (A.3) with q ú = q oe Q D .

Remark A.1

Any regular enough solution

A.3 Quadratization with the Scalar Auxiliary Variable method

Let us define a constant c oe R ú + such that 2 s L 0 U(ˆxq)dx + c is strictly positive. We also introduce an auxiliary variable which is independent of the space variable:

Associated with this auxiliary variable, we define an auxiliary operator G such that

The reformulation of the string equation can be done on the weak formulation. Seek q oe Q D et z oe R so that for all q ú oe Q D :

Since the auxiliary variable z is not distributed, the equation (A.8b) remains a non-variational scalar equation.

Any regular enough solution

Theorem A.2 (Power balance of the SAV formulation)

Proof. Use the weak formulation q ú = q oe Q D . Then equation (A.8b) allows to write that

A.4 Spatial discretization

The spatial discretization is based on high-order spectral finite elements (r oe N ú ) on a regular mesh of step x such that

The discrete finite element spaces Q h and Z h have finite dimensions n q h and n z h respectively and are defined by

and where P r is the set of polynomials of degree at most r. The degrees of freedom and quadrature points are chosen like the Gauss-Lobatto points in order to achieve mass lumping. The basis functions of Q h are denoted

, and those of Z h are denoted

. Finally, we denote Q h the vector representation of q h oe Q h and Z h that of z h oe Z h in these respective bases, and we introduce the following matrices and vectors:

where i denotes numerical integration by the Gauss-Lobatto quadrature. The auxiliary variable does not require any discretisation in space since it only depends on time. However, we give it the subscript h to indicate that it depends on the approximation in space through its coupling with the unknown q h . The semi-discretised formulation in space can then be written in matrix form:

A.5 2-SAV time discretization scheme

The variable Q h is sampled on a grid {t n = n t}, and the variable z h on an interleaved grid

inspired by the work of [Ducceschi et al., 2022], [Ducceschi and Bilbao, 2022]. The finite di erences discretisation operators are: µq n+1/2 = q n+1 + q n 2 , "q n+1/2 = q n+1 ≠ q n t , "µq n = q n+1 ≠ q n≠1 2 t " 2 q n = q n+1 ≠ 2q n + q n≠1 t 2 , {q} n ◊ = ◊q n+1 + (1 ≠ 2◊)q n + ◊q n≠1

µz n = z n+1/2 + z n≠1/2 2 , "z n = z n+1/2 ≠ z n≠1/2 t Inria