Agathe Chouippe

Eric Climent

Dominique Legendre

Céline Gabillet

Numerical simulations of drag modulation by microbubbles in a turbulent Taylor-Couette flow

Keywords: Bubbles, Taylor-Couette, Turbulence, Euler-Lagrange simulation, Two-way coupling

The aim of our study is to investigate numerically the interaction between a dispersed phase composed of microbubbles and a turbulent Taylor-Couette flow (flow within the gap between two cylinders). We use the Euler-Lagrange approach based on Direct Numerical Simulation of the continuous phase flow equations and a Lagrangian tracking for the dispersed phase. Each bubble trajectory is calculated by integrating the force balance equation accounting for buoyancy, drag, added-mass, pressure gradient, and the lift forces. The numerical method has been adapted in order to take into account the feed-back effect of the dispersed bubbles on the carrying flow. Our approach is based on local volume average of the two-phase Navier-Stokes equations. Local and temporal variations of the bubble concentration and momentum source terms are accounted for in mass and momentum balance equations. A number of reference cases have been tested to validate the modelling approach and its numerical implementation. Then, our previous study of bubble dispersion has been extended to two-way coupling simulations of turbulent Taylor-Couette flows (only inner cylinder is rotating). Modulation of the drag will be discussed for different geometries, Reynolds numbers and bubble sizes. The results show that near-wall turbulent structures are modified by the presence of bubbles.

Introduction

Studying the interaction between a dispersed phase and a continuous one has a wide variety of application, from industrial process engineering to environmental phenomena. Bubble columns, two-phase heat exchangers and bubble entrapment due to breaking waves are such examples. The aim of our study is to focus on the interaction of a dispersed phase composed of microbubbles and a turbulent flow. For that purpose we consider the configuration of the Taylor-Couette device (flow within the gap between two concentric cylinders). This flow configuration has already been investigated experimentally by [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylinder[END_REF], [START_REF] Smith | Turbulent Couette flow between concentric cylinders at large Taylor numbers[END_REF], [START_REF] Lathrop | Transition to shear-driven turbulence in Couette-Taylor flow[END_REF], and more recently at high Reynolds numbers by Van Gils et al. (2011a-c, 2012) and [START_REF] Huisman | Ultimate Turbulent Taylor-Ciyette Flow[END_REF]. Direct Numerical Simulations have been carried out for by [START_REF] Dong | Direct numerical simulation of turbulent Taylor-Couette flow[END_REF] and [START_REF] Bilson | Direct numerical simulation of turbulent Taylor-Couette flow[END_REF], and for Reynolds number up to 30000 by [START_REF] Brauckmann | Direct Numerical Simulations of local and global torque in Taylor-Couette flow up to Re=30000[END_REF]. Bubbly Taylor-Couette flows have also been experimentally explored by [START_REF] Djeridi | Bubble capture and migration in Couette-Taylor flow[END_REF] and [START_REF] Djeridi | Two-phase Couette-Taylor flow: Arrangement of the dispersed phase and effects on the flow structures[END_REF] who focused on the accumulation pattern of the dispersed phase. Dispersion mechanisms have been modelled and numerically investigated by [START_REF] Climent | Preferential accumulation of bubbles in Couette-Taylor flow patterns[END_REF] who proposed a theoretical analysis to predict the migration behaviour of bubbles and a variety of flow regimes. The feed-back modifications induced by the presence of bubbles have been characterized by [START_REF] Mehel | Analysis of the flow pattern modifications in a bubbly Couette-Taylor flow[END_REF] for moderate Reynolds numbers. Modifications of the flow structures due to bubble injection were discussed. In addition to flow structures modification, drag modulation is a key issue for many application related naval transport. For moderate Reynolds numbers, [START_REF] Murai | Frictional drag reduction in bubbly Couette-Taylor flow[END_REF] and [START_REF] Sugiyama | Microbubbly drag reduction in Taylor-Couette flow in the wavy vortex regime[END_REF] investigated experimentally and numerically the drag modification due to bubbles. This analysis has been extended to higher turbulent regimes by [START_REF] Van Den Berg | Drag reduction in Bubbly Taylor-Couette turbulence[END_REF] and Van Gils et al. (2011a-b). The aim of our study is to investigate numerically the dispersion of bubbles in turbulent Taylor-Couette flows (Reynolds number ranging from 1000 to 8000) and the possible effects on drag modulation. We are using the finite volume code Jadim with an Euler-Lagrange approach. Our paper is organized as follows: the first section is devoted to a brief description of the numerical method. We will be focusing on the model accounting for the feed-back effect of the bubbles on the carrying phase. We will also detail some test cases which have been considered to validate each physical contribution to the equations and its correct implementation in the code. The second objective of the paper is related to bubble dispersion and flow modifications. Finally, the discussion section summarizes our results in light of studies on bubbly Taylor-Couette flow, but also on bubble drag reduction observed in other geometric configurations (turbulent channel flows, turbulent boundary layers, …)

I -Numerical method

The numerical approach is based on Euler-Lagrange representation: the continuous phase flow is predicted through direct solution of the Navier-Stokes equations while Lagrangian bubble trajectories are computed by numerical integration of momentum balance equation. Modelling the presence of bubbles is based on volume-averaging of the momentum and continuity equations including momentum source terms together with spatial and temporal void fraction variations.

Equations for single phase flow

Assuming constant physical properties (viscosity and density), the Navier-Stokes equations (1.1-1.2) for a Newtonian incompressible fluid are discretized on staggered nonuniform grid with a finite volume approach.

(1.1) (1.2) Spatial derivatives are calculated with second order accuracy and we use semi-implicit Crank-Nicolson scheme for the viscous term and three steps Runge-Kutta scheme for time integration. The code has already been widely used and validated for laminar and turbulent configurations in single-phase configurations [START_REF] Climent | Preferential accumulation of bubbles in Couette-Taylor flow patterns[END_REF], [START_REF] Calmet | Statistical structure of high-Reynolds number turbulence close to the free surface of an open-channel flow[END_REF]).

Lagrangian tracking of bubbles

Trajectories are computed through temporal integration of the force balance acting on each bubble [START_REF] Climent | Large-Scale simulations of bubble-induced convection in a liquid layer[END_REF], and [START_REF] Legendre | Influence of gravity upon bubble distribution in a turbulent pipe flow: comparison between numerical simulations and experimental datas[END_REF]).

This balance takes into account buoyancy, drag, added-mass, pressure gradient, and the lift forces (eq 1.3).

v x = dt d (1.3) (1.3)
Bubbles are assumed to be spherical while the Weber number is O(10 -2 ). Drag coefficient is then estimated through the correlation of [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF], which gives the evolution of with the Reynolds number of each bubble

(1.4)

The added-mass coefficient is . And we used the correlation of [START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF] for estimating the lift coefficient variation with the bubble Reynolds number and the local shear rate. Direct bubble/bubble interactions are neglected and bubble/wall overlap is prevented by assuming an elastic bouncing. The computation of these forces requires the interpolation of the fluid velocity and its time and space derivatives at each bubble location. We used 2 nd order accuracy linear interpolation scheme. A second-order Runge-Kutta scheme is used for integrating the force balance in time with a typical time step equal to one fifth of the viscous relaxation time of the bubbles. Loop nesting is used when time steps of both phases are widely separated. When the Eulerian time step (based on numerical stability criteria) is smaller than the Lagrangian characterisitic time then solutions of both sets of equations are synchronized. For large ratio of these time steps , Lagrangian inner loops are integrated in a frozen flow field and the global computing time is dominated by the Lagrangian solver for O(10 5 to 10 6 ) bubbles. Typically, we have .

The Two-Way Coupling model

To account for the action of the dispersed phase on the carrying fluid, we average the Navier-Stokes equations in a control volume of fluid populated by bubbles. The basic principles for this operation can be found in the papers by [START_REF] Ferrante | On the physical mechanisms of drag recution in a spatially developing turbulent boundary layer laden with microbubbles[END_REF], or [START_REF] Jackson | A fluid mechanical description of fluidized beds[END_REF] and [START_REF] Drew | Mathematical modeling of two-phase flow[END_REF]. The averaging procedure is based on the function characterizing the presence of the continuous phase locally defined as (1.5) So, the presence of the dispersed phase can be defined as . We also introduce the Dirac function that enables to locate the interface within the volume such as

(1.6) (1.7)

We can also define the corresponding void fraction

(1.8)

For any given field we introduce its local phase-averaged (equation 1.9).

(1.9)

The spatial and temporal derivatives of these local-averaged function has then the following properties

(1.10)

(1.11)

We introduce the local-averaged velocity field of the continuous phase The volume averaging process leads to new continuity (1.12) and (1.13) momentum equations. The influence of the dispersed bubbles on the continuous phase is related to the local evolution of the void fraction, and to the influence of momentum source term .

(1.12)

(1.13)

(1.14)

In those equations, we have neglected momentum transport by small scale fluctuations (scales smaller than the control volume).We obtained these equations by neglecting this tensor. This contribution would be similar to the Reynolds stress tensor and may have two physical origins: turbulent subgrid scales (unresolved turbulent fluctuations), or bubble induced fluctuations due to finite size effects and direct hydrodynamic interactions. The former can be easily neglected while we intend to perform direct numerical simulations. The latter will be investigated in further studies.

Validations

We present now the different configuration tests to validate the implementation of the two-way coupling in our code. First the momentum forcing term itself has been tested without void fraction effects (case 1) , and then the modification of continuity and momentum balance equations due to non uniform (cases 2 to 6). (1.17)

We consider that drag experienced by each bubble is similar which gives:

(1.18) By integrating in the volume of the channel the force balance written in the direction we obtain an analytic expression between the pressure difference, the channel section and the sum of the forcing terms (equation 1.19).

(1.19)

The pressure has a linear evolution across the bubble cloud. While drag force is uniform, the pressure gradient must be proportional to the global void fraction as described by [START_REF] Riboux | Hydrodynamique d'un essaim de bulles en ascension[END_REF]. This is exactly what we obtained on pressure profiles (figure 3 and table 1). The agreement is very good and validates the implementation of the momentum source term. In order to validate the modification of the Navier-Stokes due to the presence of in the balance equations, we have chosen to progressively increase the complexity of tests (summarized in Table 2). For all cases, we proceed as following: we consider the configuration of case 1 (same mesh, as well) and impose simple variations of , in time or in space. So, velocity and pressure evolutions can be easily predicted by theory and compared to simulations. In cases 2 and 3 we tested the modification of the continuity equation. The aim of case 4 was to validate the modification of the pressure gradient, the modification of the advective term for case 5 and the modification of the viscous contribution to the momentum equation for case 6. For all cases we have . For cases 2 to 5, we neglected the viscous contribution while for case 6 it has been accounted for with . Figure 4 gives the evolutions of the volume fraction of the fluid, the velocity and the pressure. Table 3 summarizes the errors of the computed fields compared to analytic solutions. The error is always below 10 -10 for the velocity and O(10 -4 ) for the pressure which validates the numerical implementation.

Case

Case

Max relative error for U Max relative error for P Max error for V

Max error for W 2 3 4 5 6 Table 3 : Error obtained for the different test validation

Conclusion on the numerical method

The effect of fixed bubbles has been successfully tested. It is important to recall that the void fraction is calculated by local summation of the bubble present in a particular Eulerian mesh volume. Void fraction is transported by the Lagrangian motion of the bubbles. Finally, the code has been parallelized with distinct strategies for the Navier-Stokes solver and for the Lagrangian tracking. Domain decomposition is used for the Navier-Stokes solver thanks to [START_REF] Hallez | Mélange gravitationnel de fluids en géométrie confine[END_REF]. Bubble population is equally distributed among all processors involved in the simulation. As it will be presented in the following section, bubbly Taylor-Couette flows are characterized by local accumulation of bubbles. Distributing bubbles based on their position would lead to unbalanced processor load. Therefore, they are distributed based on their index number ranging from 1 to10 6 .

(1) x P P theo case 2 P simu case 2 P theo case 3 P simu case 3 P theo case 4 P simu case 4 P theo case 5 P simu case 5 P theo case 6 P simu case 6

(3)

0 5 10 0.1 0.15 0.2 0.25 0.3 0.35 x U x U theo x case 2 U simu x case 2 U theo x case 3 U simu x case 3 U theo x case 4 U simu x case 4 U theo x case 5 U simu x case 5 U theo x case 6 U simu x case 6

II -The Taylor-Couette flow configuration

We focus now on the simulation of turbulent Taylor-Couette flow seeded with microbubbles. We consider that only the inner cylinder is rotating and the outer one is kept fixed. The geometry (fig. 5) is characterized by the radius ratio , and the aspect ratio Γ=L ax /(R 2 -R 1 ) where L is the length of the domain, R1 and R2 are the radii of the inner and outer cylinders respectively. We use periodic boundary conditions in the axial and azimuthal directions and Dirichlet boundary conditions on cylinder walls. We will consider three radius ratios , and aspect ratios equal to for the two smallest values of η, and for the thin gap.

(2.1)

The Reynolds number (eq. 2.1) is defined with the gap width and the velocity of the inner cylinder . This dimensionless quantity which controls the dynamics of the flow for a given geometry is ranging from 1000 to 8000 in our study. For the different radius ratios and Reynolds numbers velocities are normalized by the velocity of the inner cylinder, and lengths by the radius of the outer cylinder. For even higher Reynolds number, the Taylor Vortices disappear. We will be considering the TTVF regime.

This section is organized as follows: first we describe some reference results on the single phase flow. Then, we move on the description of "passive" bubble dispersion to highlight accumulation mechanisms possibly leading to drag modulation when two-way coupling is enabled.

Single phase TTV flow

The simulations have been validated by comparing statistical quantities to previous studies. The torques acting on cylinders were compared to experimental correlations proposed by [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylinder[END_REF]. It is defined by equations

(2.1-2.2) (2.2) (2.3)
Our results of the torque are in good agreement with [START_REF] Wendt | Turbulente Strömungen zwischen zwei rotierenden konaxialen Zylinder[END_REF]. More precisely, we simulated the configuration investigated by [START_REF] Dong | Direct numerical simulation of turbulent Taylor-Couette flow[END_REF], corresponding to . Radial evolution of the mean azimuthal velocity and of the azimuthal velocity fluctuations were also in good agreement with Dong results. Averages are performed over the axial and azimuthal directions and also over time. The fluctuation is defined by . We observed that the mean velocity profiles were not very sensitive to the spatial discretisation while the fluctuation intensity was used to achieve mesh convergence. The results of our simulations are presented in fig. 7. We observe that a coarse grid leads to an overestimate of the maximum of RMS values. We also compared our results with experimental investigations of [START_REF] Mehel | Etude expérimentale d'un écoulement diphasique de Taylor-Couette[END_REF] for the case . Based on those results we selected appropriate meshes for each case for the larger gap, for the intermediate gap and for the smaller one. In order to characterize the Taylor-Couette flow it is also useful to introduce another average procedure yielding which is only averaged in the azimuthal direction and in time, as proposed by [START_REF] Bilson | Direct numerical simulation of turbulent Taylor-Couette flow[END_REF]. This allows separating the presence of large scale structures (reminiscent of Taylor vortices) and small scale turbulent fluctuations (fig. 6). The decomposition shows the specific contribution of Taylor-Vortices to RMS fluctuation (fig. 7). The presence of these Taylor Vortices induces jet zones: the Inflow and Outflow directed towards and outwards the inner cylinder. Taylor vortices and small scale fluctuations have the same intensity (between to ) for η=0.5. For thinner gaps, the ratio of the Taylor vortices contribution to the small scales contribution is obviously increased .

(1) Looking more closely at near-wall turbulent structures, we can see the presence of herringbone-like streaks. They are more intense in the outflow jet zone. These streaks have been numerically studied by [START_REF] Dong | Herringbone streaks in Taylor-Couette turbulence[END_REF], and their structure imposes the local distribution of wall shear stresses (fig. 8). Moreover, an increase of the shear stress is clearly visible in the region of jet impact (either on the inner or outer cylinder).

0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 (r -R 1 )/(R 2 -R 1 ) < U θ > xθt(r) η = 0.5 η = 0.72 η = 0.91 (2) 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 (r -R 1 )/(R 2 -R 1 ) < (U ′ θ ) 2 > xθt η = 0.5 η = 0.72 η = 0.91 (3) 0 0.2 0.4 0.6 0.8 1 0 0.05 0.1 0.15 0.2 (r -R 1 )/(R 2 -R 1 ) < (U ′′ θ ) 2 > xθt η = 0.5 η = 0.72 η = 0.91
(1) We also computed the total shear stress and tested the conservation of the kinetic momentum as: ~ as described by [START_REF] Eckhardt | Torque scaling in turbulent Taylor-Couette flow between independently rotating cylinders[END_REF] .

(2.4) Two main contributions are present: turbulent friction which dominates in the core region and viscous friction which dominates in the near wall region (fig. 9).

(1) (2)

0 0.2 0.4 0.6 0.8 1 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 (r -R 1 )/(R 2 -R 1 ) Frottement Re =
0 0.2 0.4 0.6 0.8 1 -1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 (r -R 1 )/(R 2 -R 1 ) Frottement Re =

Bubble dispersion

Bubbly Taylor Couette flow has been numerically investigated by [START_REF] Climent | Preferential accumulation of bubbles in Couette-Taylor flow patterns[END_REF] for passive dispersion.

They observed accumulation in the Taylor Vortices, which is similar to the experimental evidences of [START_REF] Djeridi | Bubble capture and migration in Couette-Taylor flow[END_REF]. Two dimensionless parameters were defined to characterize the respective contributions of the different migration mechanisms (eq. 2.5 and 2.6).

(2.5)

(2.6)

The first parameter (eq. 2.5) depends on the characteristic velocity of the Taylor Vortices and on the terminal rising velocity of the bubbles. This is a measure of force balance in the axial direction: small C values imply that bubbles are not sensitive to the Taylor Vortices and then rise in the domain. On the contrary high C values imply that buoyancy effects are small and bubbles trajectories are strongly affected by the Taylor Vortices (favourable condition for bubble entrapment). H gives an estimate of the relative influence of attraction in the Taylor Vortices in the radial direction (with high H) and centripetal attraction towards the inner cylinder due to added mass effects (typically small H). It is important to recall that H depends only on the flow parameters , while C also depends on bubbles characteristics. In the literature for higher Reynolds numbers, it has been observed that millimetric bubbles are trapped in the vortices while submillimetric bubbles are trapped near the inner cylinder, in the outflow region [START_REF] Mehel | Analysis of the flow pattern modifications in a bubbly Couette-Taylor flow[END_REF]). [START_REF] Murai | Frictional drag reduction in bubbly Couette-Taylor flow[END_REF] observed different accumulation patterns: either the formation of collars, or the formation of spirals due to buoyancy effects. In our simulations the dispersed phase is composed of 10 6 spherical bubbles, and we considered two bubble radii , two Reynolds numbers and focused on the larger gap . Bubble positions are randomly seeded in a cylindrical slab located close to the inner cylinder (between and ). Bubble velocities are initialized at the local fluid velocity. We obtained bubble accumulation only for the larger bubbles. Bubbles attracted towards the inner cylinder accumulate in the herringbonelike streaks (fig. 11) when buoyancy has a minor effect. With small bubbles we did not observe any accumulation patterns. Small bubbles are actually much more influenced by small scale turbulent fluctuations and accumulate neither in the Taylor Vortices nor in the outflow region. For large bubbles, accumulation is clear on the mean void fraction profiles (fig. 12) which is peaking near the inner cylinder. Bubble concentration maxima near the inner cylinder have also been experimentally observed by [START_REF] Mehel | Analysis of the flow pattern modifications in a bubbly Couette-Taylor flow[END_REF] and [START_REF] Murai | Frictional drag reduction in bubbly Couette-Taylor flow[END_REF] and numerically by [START_REF] Sugiyama | Microbubbly drag reduction in Taylor-Couette flow in the wavy vortex regime[END_REF]. Such bubble concentration profiles point out the limit of our modelling approach overestimating peak values compared to [START_REF] Mehel | Analysis of the flow pattern modifications in a bubbly Couette-Taylor flow[END_REF] and [START_REF] Murai | Frictional drag reduction in bubbly Couette-Taylor flow[END_REF]. ).

Indeed direct bubble interactions (such as collision) may prevent strong accumulation and helps dispersing the bubble by turbulence. As a preliminary conclusion, we can say that bubbles are sensitive to the three major contributions of the Turbulent Taylor Vortex Flow: the primary azimuthal flow induce an important attraction towards the inner cylinder, the Taylor vortices force accumulation (either in vortices or outflow regions), and the small scale turbulent fluctuations enhance local mixing of the dispersed phase. Small bubbles are far more sensitive to these small scale turbulent fluctuations, which leads to an overall uniform distribution of the void fraction.

(1) (2)

It is finally important to notice that bubbly Taylor Couette flow also presents accumulation in the radial direction for high Reynolds number cases [START_REF] Van Gils | Highly Turbulent Taylor-Couette flow[END_REF]), but it is difficult to compare our simulations with the high turbulent regime since the structure of the carrier flow is different than in our simulations.

Two-way coupling simulations

We account for the feedback effect of the dispersed phase on the carrying fluid. These aspects have been experimentally investigated by [START_REF] Murai | Frictional drag reduction in bubbly Couette-Taylor flow[END_REF] for Reynolds numbers ranging from 600 to 4500. They obtained noticeable drag reduction due to the presence of bubbles for moderate Reynolds numbers. [START_REF] Sugiyama | Microbubbly drag reduction in Taylor-Couette flow in the wavy vortex regime[END_REF] numerically investigated the bubbly drag reduction for Reynolds numbers ranging from 600 to 2500. They observed drag reduction when coherent structures of the flow were disorganized by the bubbles. Modifications of the Taylor vortices by the dispersed phase (increase of the axial wavelength of the Taylor cells) have also been experimentally observed by [START_REF] Mehel | Analysis of the flow pattern modifications in a bubbly Couette-Taylor flow[END_REF]. Moreover, bubble drag reduction has been experimentally investigated at high Reynolds numbers by [START_REF] Van Den Berg | Drag reduction in Bubbly Taylor-Couette turbulence[END_REF] and more recently in the T 3 C device by Van Gils et al. (2011).

For such high Reynolds numbers, deformability of bubbles can play a significant role. Drag reduction has also been investigated in other geometrical configuration, such as [START_REF] Ferrante | On the physical mechanisms of drag recution in a spatially developing turbulent boundary layer laden with microbubbles[END_REF] for a spatially developing turbulent boundary layer. Using a modelling approach very much similar to ours, they emphasized role of effective compressibility: bubble migration near the wall induce a non-zero divergence of the flow field which moves the turbulent structures away from the wall, leading to a global friction reduction.

In our study, Weber numbers have been evaluated with the data of the passive bubble dispersion and are typically equal to 10 -2 . Spherical shape of the bubbles can be safely assumed.

To study two-way coupling, we chose to consider the case with the smallest bubbles. Actually, simulations with large bubbles produced strong numerical instabilities due to strong local accumulation of bubbles in near-wall turbulence structures (a model to prevent bubble overlapping is required for such conditions). The dispersed phase is composed of 10 6 bubbles, leading to an average void fraction . Simulations with these characteristics did not cause any modification of the carrier flow. By multiplying the forcing terms by a factor we can reach an average void fraction . Each bubble represents a cluster of 1000 bubbles. This can be acceptable if the number of simulated bubbles is sufficient to statistically describe the dynamics of the dispersed phase. With this configuration, we obtain some modifications of the structure of the flow (fig. 13) that are similar to observations of [START_REF] Sugiyama | Microbubbly drag reduction in Taylor-Couette flow in the wavy vortex regime[END_REF] : enhancement of small scale structures and reduce strength of the Taylor vortices. We did not observe significant modification of the mean velocity profiles. We observe two different transient trends (fig. 14): drag increase on the inner cylinder and drag reduction on the outer one. Such modulation is related to the evolution of the turbulent shear stress (fig. 15) which presents a decrease in the cross-correlation of the fluctuating velocity of the coherent motion as described by [START_REF] Sugiyama | Microbubbly drag reduction in Taylor-Couette flow in the wavy vortex regime[END_REF]. Simulations must be continued to get a definite answer on whether bubbles are reducing or increasing drag in this range of parameters. Nevertheless, similar results to the study of [START_REF] Sugiyama | Microbubbly drag reduction in Taylor-Couette flow in the wavy vortex regime[END_REF] confirm that the major role of bubble modulation is located in the core of the flow.

A possible interaction scenario would be an interaction between bubbles and the inflow-outflow jets. This is still hypothetical at this stage of the work. The study of the single-phase case highlighted that high friction zones correspond to impacting jets. The rotation of the inner cylinder induces an important centripetal attraction towards the inner cylinder. Bubble migration may increase the impact of the inflow jet and decrease the intensity of the outflow, leading to an increase of the inner torque and reduction of the outer torque. Our simulations for single-phase flows showed that the relative influence of the jets tends to decrease with the Reynolds number. [START_REF] Murai | Frictional drag reduction in bubbly Couette-Taylor flow[END_REF] observed crossover between drag reduction and drag increase for Reynolds numbers near from . It would be then interesting to consider other Reynolds numbers, in order to investigate more systematically those interactions.

III -Conclusions and future prospects

We have implemented a numerical approach which enables to compute bubble dispersion with two-way coupling based on Euler-Lagrange approach. The formulation of the influence of the dispersed phase on the continuous one is based on local averaging that modifies the Navier-Stokes equations. The main modifications are the presence of a forcing term in the momentum balance equation and effective compressibility effects. Our study consists in four contributions: the first one concerns the validation of the feed-back effect on simple configurations. We obtained a very good agreement between computed and analytic solutions. The second part of our study is related to the Taylor-Couette single phase flow. Turbulent Taylor-Couette flow, in the range of considered Reynolds numbers, is composed of three major contributions: a primary azimuthal flow induced by the rotation of the inner cylinder, large scale Taylor Vortices inducing two jet zones, and small scale turbulent fluctuations. The third part of our study deals with the analysis of passive bubble dispersion. We observed that bubble migration is controlled by the three contributions of the flow. The azimuthal flow induces a sustained migration towards the inner cylinder, the Taylor Vortices may generate accumulation in the axial direction and turbulence damps local accumulation. The last part of our study concerns turbulence modulation induced by bubbles. We noticed a global decrease of turbulent fluctuations located in the core of the flow, a rearrangement of the flow in the near wall region while the mean flow profile remains unchanged. Concerning the drag modulation, we observed different trends for the inner and outer cylinders. In order to conclude more precisely on reduction or increase of drag, it will be necessary to consider other configurations, mainly by changing the global void fraction and the Reynolds number of the flow. Moreover, statistics on two-way coupling simulations need to be further converged.

Figure 1 :

 1 Figure 1 : Schematic representation of a fluid control volume populated with bubbles, corresponds to the interface between both phases We introduce then the local volume fraction of the continuous phase (equation 1.7).

  In this first test-case we consider a 3D-channel with stress free boundary conditions on the four boundaries parallel to the mean flow (symmetry conditions) and inflow-outflow conditions in the axial direction. In this case bubbles remain fixed at random positions in the central part of the simulation domain (fig.2). The dimensions of the channel are and . The mesh is regular in x and z directions and non-uniform in y direction with , We impose a constant inlet velocity with . Simulations were unsteady. Pressure is obtained by an iterative solver for the Poisson equation with a convergence criterion of . Fluid characteristics are chosen to lead to negligible viscous dissipation.

Figure 2 :

 2 Figure 2 : Sketch of the configuration used for the validation in the case 1 whith 40,000 uniform bubbles of radius

Figure 3 :

 3 Figure 3 : Evolution in the direction of the pressure averaged in the channel section, for 3 different void fractions (a) and scaled by the global void fraction (b)

  Characteristics of the validation tests for non-uniform void fraction distribution.

Figure 4 :

 4 Figure 4 : Summary of the fluid concentration (1), pressure (2) and velocity (3) for the different cases over in the x direction. For case 3 it corresponds to time

Figure 5 :

 5 Figure 5 : Sketch of the Taylor-Couette configuration.For low Reynolds numbers the flow is purely azimuthal. When increasing we have the TVF regime (Taylor Vortex flow) characterized by the development of large scale toroidal counter-rotating vortices. For larger Reynolds number the Wavy Vortex Flow (WVF) regime occurs followed by Modulated Wavy Vortex Flow (MWVF). At larger Reynolds numbers, the flow becomes turbulent (Turbulent Taylor Vortex Flow, (TTVF), which is characterized by the persistence of the large scale Taylor Vortices and the presence of turbulent fluctuations (fig.6). For even higher Reynolds number, the Taylor Vortices disappear. We will be considering the TTVF regime.

Figure 6 :

 6 Figure 6 : Isosurfaces of the velocity field for , (Left: Instantaneous flow field, right: Mean velocity field )

Figure 7 :

 7 Figure 7 : Statistics on the azimuthal component on single-phase flow for and for different radius ratios: (1) mean field, (2) RMS fluctuations , (3) RMS fluctuations (((u" θ ) 2 ) θt ) 1/2 .

Figure 8 :

 8 Figure 8 : Wall shear stresses for and . (1) Inner cylinder, (2) Outer cylinder

Figure 9 :

 9 Figure 9 : Different contributions of the global shear for and for different Reynolds numbers: (1) turbulent shear, (2) viscous shear The pdf of the local friction (fig. 10) is not symmetric (similarly to turbulent channel flows).(1)

Figure 10 :

 10 Figure 10 : Normalized pdf of the wall shear stress for and for different Reynolds numbers: (1) for the inner cylinder, (2) for the outer cylinder

Figure 11 :

 11 Figure 11 : Bubble positions near the inner cylinder for the case , , H=0.0637 and Top: without buoyancy, bottom: with buoyancy ().

  Figure 12 : Mean void fraction profiles at different times for case , , H=0.0637, C=∞: (1) for(2)

Figure 13 :

 13 Figure 13 : Isovalue of the velocity field ( ) after three rotations of the inner cylinder (colors corresponds to the axial velocity, and highlights the presence of the Taylor Vortices). Left: one-way coupling, right: two-way coupling.

Figure 14 :Figure 15 :

 1415 Figure 14 : Temporal evolution of the torques acting on the inner cylinder (1) and outer cylinder (2). Time is scaled with the one revolution time of the inner cylinder

Table 1 :

 1 Pressure gradient and relative errors obtained for the different cases considered for the validation of the implementation of the momentum forcing term

b)

Cases 2 to 6: modifications of the Navier-Stokes equations by a non-uniform void fraction

Table 4

 4

	summarizes the parameters we used. We focus on cases
	without gravity. Wall units are based on the friction velocity
	for each wall		, the associated viscous length
	scale	and turbulent Reynolds number
		.
			η 0.5
		Re 5000
		/
		( -)/	1250
			0.23
			C ∞
			H 0.0637
			300
			150

Table 4 :

 4 Parameters of the two-way coupling simulations.

Acknowledgements

We would like to acknowledge Délégation Générale de l'Armement (DGA) for funding this project. The regional supercomputing centre CALMIP at the University of Toulouse, and the national centres IDRIS and CINES are gratefully acknowledged for their support. We also thank Annaïg Pedrono for her technical support.


C Irimea

I K Ortega

Y Carpentier

Michael Ziskind

F X Ouf

F Salm

D Delhaye

D Gaffié

A Faccinetto

C Focsa

Airplane

Airplane Soot Chemical Composition Analysis by Laser Desorption and Secondary Ions Mass Spectrometry

Keywords: soot, PAH, aviation, laser desorption, SIMS

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Soot emissions from airplane engines affect both air quality and climate. Up to 10% increase of cirrus clouds formation has been observed in flight corridors during the last decade (Boucher 1999). So most likely soot particles emitted by air traffic allow the formation of ice at lower supersaturations compared to homogeneous ice formation. The potential of soot particles as cloud condensation nuclei is affected by their physico-chemical properties, for example the presence of Polycyclic Aromatic Hydrocarbons (PAHs) adsorbed on the soot surface plays a role in its reactivity and ice nucleation potential (Crawford et al. 2011).

The present work has been performed in the framework of the French national project MERMOSE, which aims to characterize the soot emissions from airplane engines and to determine their role in the formation of cirrus clouds. A first sampling campaign performed on the PowerJet SaM-146 turbofan (SNECMA/Saturn) provided soot particles collected on quartz fiber filters in four engine regimes (30%, 70%, 85% and 100%). Two-Step (Desorption/Ionization) Laser Mass Spectrometry (L2MS) and Time of Fight Secondary Ion Mass Spectrometry (ToF-SIMS) techniques were used to analyze the surface chemical composition of the soot particles, with special focus on PAHs. In L2MS, the adsorbed phase is probed by nanosecond laser desorption (λ d =532nm), then the ejected molecules are ionized with a second ns laser (λ i =266nm) and further mass-separated by ToF-MS.

Figure 1. Positive polarity mass spectrum of a soot sample collected in the 100% engine regime.

Excellent sensitivity (limit of detection below 1 fmol per laser pulse) is reached by this technique especially for PAHs, thanks to their large REMPI (Resonance-Enhanced Multi-Photon Ionization) crosssections at 266 nm. Figure 1 shows an example of a typical mass spectrum obtained during the L2MS experiments. In ToF-SIMS the sample is bombarded with a Bi 3 + ion beam and the secondary ions generated are detected by ToF-MS. In both techniques the spectra were obtained using positive polarity, which is better suited for detection of PAHs. A good agreement was obtained between the two techniques for the total PAH content of the analyzed samples.

Moreover, the total PAH content followed the same trend as the OC/EC ratio measured with a thermo-optic analyzer (Improve protocol): the 30% engine regime soot presents a high concentration of PAHs and a high OC content, while the three other regimes give a relatively low content of PAHs and OC. Figure 2 summarizes these results. surrogate lab generated soot. Furthermore, the ice nucleation properties of these samples will be measured using Raman spectroscopy.

Figure 2 .

 2 Figure 2. Left panel: total PAHs content of different samples analyzed by TOF-SIMS (dark green) and L2MS (light green) techniques; right panel: organic carbon (dark green) and elemental carbon (black) content of different samples analyzed.This study is being extended to other compounds (e.g. organosulfates or oxygenated hydrocarbons) and to

Acknowledgments

This work was supported by the French Civil Aviation Authority (DGAC) through the MERMOSE project (http://mermose.onera.fr/).


Emmanuel Giguet

Nadine Lucas

come


Juliette Chabassier
email: juliette.chabassier@inria.fr

Marc Duruflé
email: marc.durufle@inria.fr

Energy based simulation of a Timoshenko beam in non-forced rotation. Influence of the piano hammer shank flexibility on the sound

Keywords: flexible piano hammer shank, piano model, pianistic touch, implicit time discretization, vibrating Timoshenko beam, energy preservation 2011 MSC: 35L05, 35A35, 35A40, 74K05

Energy based simulation of a Timoshenko beam in non-forced rotation. Application to the flexible piano hammer shank.

Motivation

Experimental studies [START_REF] Koornhof | The influence of touch on piano sound[END_REF][START_REF] Goebl | Once again: The perception of piano touch and tone. Can touch audibly change piano sound independently of intensity ?[END_REF] tend to conclude that "the difference between two equally loud piano tones due to type of touch lies in the different noise components involved in the keystroke" (from [START_REF] Goebl | Once again: The perception of piano touch and tone. Can touch audibly change piano sound independently of intensity ?[END_REF]). This "noise" is supposed to come from the shock excitation of the piano structure (keybed, rim, iron plate). The authors seem to exclude any other mechanism explaining an influence of the pianistic touch on the resulted piano sound. However, experimental observations of [START_REF] Askenfelt | From touch to string vibrations. II: The motion of the key and hammer[END_REF][START_REF] Giordano | Piano hammers and their force compression characteristics: Does a power law make sense ?[END_REF][START_REF] Birkett | Experimental investigation of the piano hammer-string interaction[END_REF] suggest that the vibration of the hammer shank plays a significant role in the interaction between the piano hammer and the string without being able to quantify its effect on the sound. Different possible mechanisms have been suggested, as a change in the striking position or a rubbing motion in the elongation direction of the string [START_REF] Askenfelt | From touch to string vibrations. II: The motion of the key and hammer[END_REF]. But the most frequent suggestion is that the shank vibrations contribute to the interacting hammer headstring force, and this is the mechanism that we wish to investigate in this article. Recently, numerical experiments have highlighted the influence of the hammer shank's flexibility on the whole piano action [START_REF] Izadbakhsh | Dynamic Modeling and Experimental Testing of a Piano Action Mechanism With a Flexible Hammer Shank[END_REF] and on the interaction between the hammer head and the strings [START_REF] Vyasarayani | Modeling the dynamics of a compliant piano action mechanism impacting an elastic stiff string[END_REF]. To our knowledge, there is no published work which shows numerical experiments of a flexible hammer shank coupled with a whole piano (strings, soundboard and sound radiation in the air), therefore showing the influence of the shank flexibility on the resulting sound.

The aim of this work is to provide a complete physics based modeling and numerical solution of the interaction between a flexible hammer shank and a hammer head, using only constitutive laws. Our goal is also to design an energy-based numerical solution method, such that the numerical methods will ensure the preservation (or decay) of a discrete energy approximating the mechanical energy of the continuous system. Coupling this new numerical method to the energy-preserving numerical solution for the other parts of the piano described in [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF] (beginning with the strings, but also a soundboard and sound radiation in the free space) will thus be easily possible in a stable and physically relevant way. As a result, we will be able to provide several piano sounds obtained with the exact same configuration, where only the input forces (associated to different pianistic touches) will be changed. Moreover, our simulations do not model the shock on the structure. The impact point will be exactly the same is all simulations, and none of them implement a "longitudinal rubbing motion" of the hammer head on the string. Differentiating physical phenomena is a powerful feature of numerical experiments that cannot be done in a real experimental setup. Therefore the observed differences on these sounds will be attributed to the touch only.

A rotating Timoshenko beam will represent the flexible hammer shank. The beam will be in interaction with a mass located at its extremity representing the hammer head. The "macroscopic" motion of the beam is a rotation, the rotation angle being an unknown of the problem. The "microscopic" motions of the beam are small vibrations that can be viewed as small perturbations around the macroscopic position. The interacting forces acting between the hammer head and the hammer shank will be considered as Lagrange multipliers of the system. The obtained system of equations will involve nonlinear terms that will arise from the rotation of the coordinates system. We will see that some of them can be neglected, but the problem will remain challenging to solve numerically. The mechanical system that we consider is a small part of the grand piano action (see figure 1) and consists of two parts:

• the hammer shank, in which we consider linear beam vibrations around a macroscopic position, written in rotating coordinates. It is subject to gravity, a "pushing" force (coming in reality from the contact with the jack), and a force at its extremity coming from the hammer head.

• the hammer head, not flexible, attached to the free end of the beam (cinematic continuity). It is subject to gravity, acts back on the hammer shank, and interacts with the strings.

Though the final goal of this work is to model the coupling of this mechanical system with piano strings, in the two first parts of this article we will simplify the problem and suppose that the forces coming from the coupling with strings are given data of the problem. In the third part of this article, we will provide numerical examples where the hammer shank and head have been coupled to the strings, the soundboard and the sound radiation, by treating the hammer head / string coupling in a mathematical weak way using Lagrange multipliers and computing the interaction forces from a crushing felt law as in [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF].

A model for the { shank -head } mechanical system

We consider the notations of figure 2(a). The origin O is located at coordinates (x 0 , y 0 , z 0 ) = (0, 0, 0) in the cartesian reference frame (u x , u y , u z ) This point will be fixed, since the only macroscopic motion of the shank that we consider is a rotation, with an angle θ.

In the following of this article, we will consider that all exterior forces are in the (u x , u y ) plane, therefore the displacement field will occur in this plane. Thus, we will only seek mechanical solutions in the (u x , u y ) plane. We denote the three unit vectors :

u r =           cos θ sin θ 0           , u z =           0 0 1           , u θ =           sin θ -cos θ 0          
The system (u θ , u r ) is orthonormal, we denote by (s, τ) the coordinates of the undeformed configuration in this system (see figure 2(b)).

Following Timoshenko theory [START_REF] Timoshenko | Theory of plates and shells, Second Edition[END_REF], which states that straight cross-sections remain straight during the deformation, the displacement of the beam's cross-section defined by {s u r + τ u θ , τ ∈ (-R, R)} is chosen of the form:

r = (s -τ sin ϕ) u r + (w + τ cos ϕ) u θ (1)
with w the deflection of the neutral fiber of the beam and ϕ the rotational angle of the cross section. The axis of rotation of the beam is equal to Ω = u z . The vectors u r and u θ are not constant over time (in the sequel we shall denote ẋ the time derivative of any quantity x):

ur = θ Ω × u r = -θ u θ , uθ = θ Ω × u θ = θ u r
And the second order time derivatives of these vectors read:

ür = -θ2 u r -θ u θ , üθ = -θ u r -θ2 u θ
The wooden part of the head is assumed to be rigid, hence the center of gravity of the mass (hammer head) at the end of the beam is parametrized by:

ξ = L u r + (w(L) -H) u θ (2)
The top of the hammer head which can be in contact with the strings is parametrized by:

ζ = L u r + w(L) -H ′ u θ
where

H ′ = H + d 0 (see figure 2(a)).
In the sequel, we shall derive the equations of motion for the system { shank, head } using the virtual work principle. This will yield a variational formulation of the problem. Then, an energy identity will be obtained.

In the following, we denote with the symbol the integration over [0, L] and with the integration on the 3D domain occupied by the beam.

The virtual work principle states that, P acc + P int = P ext where P ext = P grav + P jack + P coupl,shank + P coupl,head + P string where P acc is the virtual work due to acceleration, P int the virtual work due to the interior forces of the shank and P ext the virtual work due to the exterior forces. This last item can be divided into several contributions coming from the gravity, the jack force, the coupling between the shank and the head, and the string force. The different terms are derived successively in Appendix A. Substituting them with their expression then yields a variational formulation of the problem, where A is the area of the beam section, and I is the bending moment of inertia of the beam. For instance, in the case of a circle cross-section, we have

A = π R 2 and I = π R 4 4
, where R is the radius.

Variational formulation of the problem

The internal coupling force F coupl is an unknown of the problem which has been introduced in order to treat the shank-head coupling in a mathematical weak way. The additional equation comes from imposing that the displacement of the hammer head is equal to the displacement of the extremity of the shank. This equation is called "equation of continuity shank-head".

The variational formulation directly yields from replacing the virtual works by their expression. In order to prepare for the numerical resolution, it is moreover necessary to specify the variational spaces of the solutions. Let us introduce

V 0 = {w ∈ H 1 ([0, L]), such that w(0) = 0}.
One major hypothesis for writing the displacement field of the shank as equation [START_REF] Koornhof | The influence of touch on piano sound[END_REF] was that all exterior forces were in the (u x , u y ) plane. Using this hypothesis, it is straightforward to show that the solution lives in the same plane. In the following, we will only consider a 2D motion (meaning that ξ and F coupl ∈ R 2 ), but for convenience we will still call u x , u y , u r and u θ their projection on the 2D plane :

u x = 1 0 , u y = 0 1 , u r = cos θ sin θ , u θ = sin θ -cos θ .
The reader can acquire a first intuition of the fact that an additional constraint on w is needed in order to grant the uniqueness of solutions, by realizing that several couples (θ, w) can describe a given beam position. For instance, (θ = 45 • , w = 0) describes the same situation as (θ = 0, w(s) = -s), as an infinity of other couples. But, there is a unique couple which satisfies that the mean value of w is zero. Of course this criterion for the choice of a unique couple is coherent with the physical assumption for using Timoshenko theory of "small deformations". An analytical study of this criterion is beyond the scope of this article and shall be addressed in an other context. We chose to impose that the mean value of w is zero via a Lagrange multiplier technique. The Lagrange multiplier is called λ.

                                                                                                                                                     We seek (w, ϕ, θ, λ, F coupl , ξ) ∈ V 0 × H 1 × R × R × R 2 × R 2
such that for all (w ⋆ , ϕ ⋆ ) ∈ H 1 × H 1 , we have Shank equations / equation on (w, ϕ):

(3a)

L 0 ρA ∂ 2 w ∂t 2 w ⋆ + L 0 AGκ ∂w ∂s -ϕ ∂w ⋆ ∂s + λ L 0 w ⋆ -ρA w θ2 + θ s + g cos θ w ⋆ -F coupl • u θ w ⋆ (L) = -F jack L 0 δ jack (s) w ⋆ (s) ds L 0 ρI ∂ 2 ϕ ∂t 2 ϕ ⋆ - L 0 AGκ ∂w ∂s -ϕ ϕ ⋆ + L 0 EI ∂ϕ ∂s ∂ϕ ⋆ ∂s -ρI θ L 0 ϕ ⋆ = 0
Zero average of w / equation on λ:

(3b) L 0 w = 0 Equation on θ: (3c) θ L 0 ρA (w 2 + s 2 ) + 2 θ L 0 ρA w ẇ - L 0 ρA s ẅ + L 0 ρA g (s cos θ + w sin θ) + L 0 ρI θ -φ = F jack L 0 sδ jack + (w(L) -H) F coupl • u r -L F coupl • u θ Head equation / equation on ξ: (3d) m H ξ = -m H g u y -F coupl -F string u y
Equation of continuity shank -head / equation on F coupl : (3e)

ξ = L u r + (w(L) -H) u θ

Energy identity

Taking as test functions (w ⋆ , ϕ ⋆ ) the time derivative of the unknowns (w, ϕ) in (3a), multiplying equation (3c) with θ and multiplying equation (3d) with ξ enables us to show the following energy identity:

dE dt = F jack θ sδ jack (s) ds - ẇ(s)δ jack (s) ds -F string ξ • u y
where the energy reads:

E(t) = 1 2 L 0 ρAw 2 ( θ) 2 + 1 2 L 0 ρA( ẇ -s θ) 2 + 1 2 L 0 ρI( θ -φ) 2 + m H 2 ξ 2 + 1 2 L 0 AGκ ∂w ∂s -ϕ 2 + 1 2 L 0 EI ∂ϕ ∂s 2 + L 0 ρA g (s sin θ -w cos θ) + L 2 ρ A g 2 + m H g (ξ • u y ) + C g
Especially, this energy is preserved if exterior forces F jack and F string are equal to zero.

Notice that it is obtained by using the time derivative of (3e). Indeed, it allowed us to simplify the terms involving F coupl , which are nothing but internal forces to the global mechanical system.

Most terms of this energy are trivially positive. Only the last terms, associated to gravity, need a detailed explanation:

                 L 0 ρA g s sin(θ) + L 2 ρAg/2 ≥ 0 (4a) L 0 ρA g w cos(θ) = 0 thanks to (3b) (4b)
The last term of the energy, m H g ξ • u y + C g , is not necessarily positive, but it would be with C g = -√ H 2 + L 2 if the shank was not flexible (which is the macroscopic motion of the beam). Since the beam vibrations are supposed to be small, we expect that there exists a constant C g close to its "macroscopic" equivalent that makes the contribution positive. This qualitative argument should be shown in a further analytic study which is beyond the scope of this article.

3. An energy-based numerical scheme for solving the { shank -head } system We are guided by the desire to preserve a discrete energy, which will provide a priori estimates on the discrete solution. The energetic approach is a powerful way to later introduce couplings with other mechanical systems, in a stable and preserving way, which is our final goal in the context of the piano, beginning with the strings.

We rely on the variational form given by the set of equations (3) and perform a conformal Galerkin method in space with high order finite elements 1 . Let us denote V h,0 ⊂ V 0 and V h ⊂ H 1 the approximate spaces given as

[0, L] = k [s k , s k+1 ] V h = u ∈ H 1 ([0, L]) such that u| [s k ,s k+1 ] ∈ P r V h,0 = u ∈ V 0 such that u| [s k ,s k+1 ] ∈ P r
where P r is the space of polynomials of degree lower or equal to r. Time discretisation is the main issue since it is absolutely not trivial to ensure numerical stability in such a nonlinear and coupled context, especially in view of the coupling with the rest of the instrument. Most of the terms are treated in a classical implicit way with "θ-schemes", the preserving class of Newmark schemes (see chapter XX of [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF]). The weighting coefficient is chosen equal to 1/4 in order to avoid any stability condition while maximizing the accuracy. A non classical scheme is proposed in order to discretize the inertial term. These choices lead to a stable scheme, and the over cost of implicit schemes is anyhow expected in such a nonlinear context.

We denote

{X} n 1/4 = X n+1 + 2X n + X n-1 4 (5)
which is a centered weighted approximation of X(t n ). In the same way, we denote

{cos θ} n 1/4 := cos θ n+1 + 2 cos θ n + cos θ n-1 4 (6)
which is a centered weighted approximation of cos θ(t n ) where the average is done outside the cosine function. We also introduce the functions

δ cos(a, b) =            cos a -cos b a -b if a b -sin(a) if a = b and δ sin(a, b) =            sin a -sin b a -b if a b cos(a) if a = b (7)
and finally, the oriented vectors

u n r = cos θ n sin θ n , u n θ = sin θ n -cos θ n
In the sequel, numerical integration is done with quadrature techniques (see [START_REF] Quarteroni | Méthodes Numériques[END_REF]) and denoted .

Non classical treatment of the inertial term

A non-classical discretisation is done for the three contributions coming from the first term of the energy:

1 2 L 0 ρA w 2 ( θ) 2 (8)
which yields the three following terms in the equations:

                 - L 0 ρA w θ2 w ⋆ for the shank's flexural equation, θ L 0 ρA w 2 + 2 θ L 0 ρA w ẇ for the equation on θ. (9)
A first idea would be to use centered weighted approximations for w and centered quotients for the time derivatives of θ and w, but this does not lead to a conservative numerical scheme. A second idea would be to use variational time integrators methods as described in [START_REF] Lew | Variational time integrators[END_REF], but these methods do not preserve the discrete energy if the time step is constant. It is however a good starting point to find energy preserving schemes. Let us call L h the discrete density of lagrangian associated to the continuous Lagrangian (8):

L h (q n , q n+1 ) = 1 2 θ n+1 -θ n ∆t 2 w n+1 + w n 2 2
, where q n = (w n , θ n ) [START_REF] Dautray | Mathematical analysis and numerical methods for science and technology[END_REF] Then the variational integration technique yields the following scheme:

                       1 2 θ n+1 -θ n ∆t 2 w n+1 + w n 2 + 1 2 θ n -θ n-1 ∆t 2 w n + w n-1 2 ← w θ2 (11a) 1 ∆t w n+1 + w n 2 2 θ n+1 -θ n ∆t - 1 ∆t w n + w n-1 2 2 θ n -θ n-1 ∆t ← θ w 2 + 2 θw ẇ (11b)
This method does not preserve exactly the discrete Lagrangian L 0 ρA L h , but we use the second line of this scheme as a starting point and change the first line in order to guarantee an energy preservation. Finally, the three contributions will be discretized as follows:

                       - L 0 ρA w n+1 + 2w n + w n-1 4 θ n+1 -θ n ∆t θ n -θ n-1 ∆t w ⋆ ← - L 0 ρA w θ2 w ⋆ (12a) L 0 ρA ∆t        θ n+1 -θ n ∆t w n+1 + w n 2 2 - θ n -θ n-1 ∆t w n + w n-1 2 2        ← θ L 0 ρA w 2 + 2 θ L 0 ρA w ẇ (12b)
The consistency of system [START_REF] Chapelle | An optimal low-order locking-free finite element method for Reissner-Mindlin plates[END_REF] with the inertial terms of the equation is not straightforward and is shown in Appendix B. Since this scheme is not classical, we provide the calculations that lead to the discrete energy.

We take w ⋆ = w n+1w n-1 2∆t in (12a) and multiply (12b) by

θ n+1 -θ n-1 2∆t
, and we add the contributions. We get, after

denoting w n+ 1 2 = w n+1 + w n 2 and θn+ 1 2 = θ n+1 -θ n ∆t , 0 = L 0 - ρA 2∆t θn+ 1 2 θn-1 2 w n+ 1 2 2 -w n-1 2 2 + L 0 ρA 2∆t w n+ 1 2 2 θn+ 1 2 2 + w n+ 1 2 2 θn+ 1 2 θn-1 2 -w n-1 2 2 θn-1 2 2 -w n-1 2 2 θn+ 1 2 θn-1 2
The two terms coming from (12a) both cancel with terms coming from (12b), and we get :

E n+ 1 2 NL -E n-1 2 NL ∆t = 0, where E n+ 1 2 NL = L 0 ρA 2 θ n+1 -θ n ∆t 2 w n+1 + w n 2 2 8

Numerical scheme

The numerical scheme that we propose, based on our innovative treatment of the inertial terms, is the following:

                                                                                                                                                                                                                   We seek (w n+1 , ϕ n+1 , θ n+1 , λ, F coupl , ξ n+1 ) ∈ V h,0 × V h × R × R 2 × R 2 such that for all (w ⋆ , ϕ ⋆ ) ∈ V h,0 × V h , we have Shank equations: (13a) L 0 ρA w n+1 -2w n + w n-1 ∆t 2 w ⋆ + L 0 AGκ ∂w ∂s -ϕ n 1/4 ∂w ⋆ ∂s + λ L 0 w ⋆ -ρA w n+1 + 2w n + w n-1 4 θ n+1 -θ n ∆t θ n -θ n-1 ∆t + θ n+1 -2θ n + θ n-1 ∆t 2 s + g {cos θ} n 1/4 w ⋆ -F coupl • u n θ w ⋆ (L) = -F jack L 0 δ jack (s)w ⋆ (s) ds L 0 ρI ϕ n+1 -2ϕ n + ϕ n-1 ∆t 2 ϕ ⋆ - L 0 AGκ ∂w ∂s -ϕ n 1/4 ϕ ⋆ + L 0 EI ∂ϕ ∂s n 1/4 ∂ϕ ⋆ ∂s -ρI θ n+1 -2θ n + θ n-1 ∆t 2 L 0 ϕ ⋆ = 0 Zero average of w: (13b) L 0 {w} n 1/4 = 0 Equation on θ: (13c)
θ n+1 -2θ n + θ n-1 ∆t 2 L 0 ρA s 2 + L 0 ρA ∆t        θ n+1 -θ n ∆t w n+1 + w n 2 2 - θ n -θ n-1 ∆t w n + w n-1 2 2        + L 0 ρA g s δ sin(θ n+1 , θ n-1 ) -{w} n 1/4 δ cos(θ n+1 , θ n-1 ) - L 0 ρA s w n+1 -2w n + w n-1 ∆t 2 + L 0 ρI θ n+1 -2θ n + θ n-1 ∆t 2 - ϕ n+1 -2ϕ n + ϕ n-1 ∆t 2 = F jack L 0 sδ jack + {w} n 1/4 (L) -H F coupl • u n r -L F coupl • u n θ Head equation: (13d) m H ξ n+1 -2ξ n + ξ n-1 ∆t 2 = -m H g u y -F coupl -F string u y
Equation of continuity shank -head: (13e)

ξ n+1 -ξ n-1 2∆t = θ n+1 -θ n-1 2∆t {w} n 1/4 (L) -H u n r + w n+1 -w n-1 2∆t (L) -L θ n+1 -θ n-1 2∆t u n θ
Notice that the equations of continuity are consistent with the time derivative of conditions (3e). This will enable us to show a discrete energy preservation.

Discrete energy

A discrete energy identity can be obtained by following the classical steps:

• taking in (13a) as test functions :

w ⋆ = w n+1 -w n-1 2∆t , ϕ ⋆ = ϕ n+1 -ϕ n-1 2∆t ,
• multiplying (13c) by θ n+1 -θ n-1 2∆t ,

• multiplying (13d) by

ξ n+1 -ξ n-1 2∆t .
Notice that by definition,

δ sin(θ n+1 , θ n-1 ) θ n+1 -θ n-1 = cos(θ n+1 ) -cos(θ n-1 ) and δ cos(θ n+1 , θ n-1 ) θ n+1 -θ n-1 = sin(θ n+1 ) -sin(θ n-1 ).
The interior coupling terms will vanish thanks to equation (13e).

Without any source term (F jack = F string = 0), the following discrete energy is preserved from one time step to the following:

E n+ 1 2 = E n+ 1 2 kin + E n+ 1 2 pot , where (14)
E

n+ 1 2 kin = 1 2 L 0 ρA θ n+1 -θ n ∆t 2 w n+1 + w n 2 2 + 1 2 L 0 ρA w n+1 -w n 2∆t -s θ n+1 -θ n 2∆t 2 + 1 2 L 0 ρI θ n+1 -θ n ∆t - ϕ n+1 -ϕ n ∆t 2 + m H 2 ξ n+1 -ξ n ∆t 2 (15)
E

n+ 1 2 pot = 1 2 L 0 AGκ ∂ ∂s w n+1 + w n 2 - ϕ n+1 + ϕ n 2 2 + 1 2 L 0 EI ∂ ∂s ϕ n+1 + ϕ n 2 2 + L 0 ρA g s sin θ n+1 + sin θ n 2 + L 2 ρAg/2 - L 0 ρA g w n+1 + w n 2 cos θ n+1 + cos θ n 2 + m H g ξ n+1 + ξ n 2 • u y + C g (16)
Using classical tools (see [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF]), the positivity of this quantity would grant the unconditional stability of the numerical scheme. As in the continuous case, we have:

                 L 0 ρA g s sin θ n+1 + sin θ n 2 + L 2 ρAg/2 ≥ 0 L 0 ρA g w n+1 + w n 2 cos θ n+1 + cos θ n 2
= 0 thanks to (13b), and assuming that

L 0 w 0 = L 0 w 1 = 0.
The only terms which are not necessarily positive is m H g ξ n+1 + ξ n 2 • u y + C g , as in the continuous case. The same qualitative argument applies, therefore this term is expected to be positive.

Numerical simulations and prospects

For the numerical simulations that follow, the hammer shank and head are coupled to one or several vibrating strings, which are coupled to a soundboard that radiates in the air, as in a real piano. Therefore, in the following, F string is not longer a given force, but comes from the interaction between the hammer head and the string(s), which is (are) now free to vibrate. At the origin of times, the piano action is at rest, meaning that θ ≈ -17 deg. The strings are at rest and occupy the horizontal line y = H + d 0 (see figure 2(a)) so that, macroscopically, the hammer head -string contact begins when the angle θ is equal to 0. This choice was motivated by several observations in real pianos, but is not universally respected in the instrument manufacture. The hammer head -string contact is treated as in [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF] which means that strong hypotheses have been made :

• the contact region on the string is a fixed interval of the string,

• the interacting force is a nonlinear function of the distance d between the hammer head and the string, with hysteresis. This nonlinear function F string is defined as (see [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF])

F string (d) = K Φ(d) + R d dt Φ(d), Φ(d) = (d 0 -d) + p ( 18
)
where p is a positive real number, d 0 is the distance between ξ and the top of the hammer head ζ at rest (see figure 2(a)), and (•) + stands for the positive part.

In the model described in reference [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF], the hammer head was only allowed a vertical motion, whereas in the present shank model an horizontal motion of the hammer head can occur. However, for the following numerical results, the distance d has been defined as the vertical distance between the hammer head and the string. This means that only the vertical position of the hammer head counts in the contact. This is a rough description of the geometry and we will discuss a possible improved model in section 4.4.

A new feature that we had to model was the let-off mechanism. In real pianos, the jack (see figure 1) pushes the shank until the jack's other end touches a fixed obstacle which makes it rock and loose contact with the shank. In a first approach, we modeled this mechanism as the following : as soon as the distance between the hammer head and the string (at rest) is smaller than the so-called let-off distance, the jack force is set to zero. In a real piano, the contact element between the jack and the shank (the roller) is covered with felt, and therefore the contact loss is probably smoother. In the absence of experimental data concerning this felt and its crushing law, we chose this abrupt solution but a better model could be implemented in the future.

The modeling of the strings, soundboard and air is beyond the scope of this article, so we refer the reader to [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF] for a description of the mechanical model of the whole piano with stiff and nonlinear vibrating strings, and to [START_REF] Chabassier | Stability and dispersion analysis of improved time discretization for simply supported prestressed Timoshenko systems. Application to the stiff piano string[END_REF] and [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF] for a description of the innovative energy preserving numerical schemes used to solve the equations. The coupling between the numerical scheme (13) and the numerical scheme done on the strings and the rest of the piano is handled via Lagrange multipliers which are the coupling forces F string (one for each string).

The numerical experiments are conducted with the shank's parameters given in table 4 (taken from [START_REF] Izadbakhsh | Dynamic Modeling and Experimental Testing of a Piano Action Mechanism With a Flexible Hammer Shank[END_REF]) and the head's parameters given in table 4. The piano's parameters can be retrieved in [START_REF] Chabassier | Physical parameters for piano modeling[END_REF]. The jack force F jack is applied around a point located 0.0155 m from the rotation center of the shank (point (0, 0) in figure 2(a)). The repartition function δ jack stands for the roller's width. The time dependency of the jack force is a data of our simulations that we need to provide. In the literature, no experimental data of grand piano jack-shank forces and hammer -string interaction forces are exploitable, partly because they seem difficult to achieve with precision and/or because of unjustified assumptions in the experimental setup (see [START_REF] Birkett | Experimental investigation of the piano hammer-string interaction[END_REF]). However, two examples of jack-shank forces can be found in [START_REF] Askenfelt | From touch to string vibrations. II: The motion of the key and hammer[END_REF] as reproduced in figure 3. From this figure, we have isolated two kinds of piano touches, trying to mimic the so-identified legato and staccato touches. The staccato touch will be associated with a strong and short pulse as a jack force, returning to zero before any let-off mechanism. The legato touch will be associated with an increasing jack force which only stops because of the presence of a let-off mechanism. Both touches will be modeled with the following formula for F jack , only changing the parameters t jack and A jack .

L (m) E (Pa) A (m 2 ) ρ (kg • m -3 ) I (m 4 ) G (Pa) κ (
F jack (t) =              A jack cos 2 π 2 2 t t jack -1 , if 0 ≤ t ≤ t jack 0, otherwise (19)
Once again, this is only a rough modeling of a very complex situation, which could be improved in the future.

The purpose of the following numerical simulations is first to numerically demonstrate the preservation of a discrete energy (see 4.1), and more interestingly, to illustrate the versatility of our simulation tool. The nonlinear system that we modeled and solved in this article depends on several physical parameters which are very interesting to explore. We will first investigate the effect of the pianistic touch via the jack force (see 4.2), and of the let-off distance (see 4.3). During these investigations, we will discuss the possible presence of a "longitudinal touch precursor" (see 4.4) in the piano sound.

Numerical algorithm and energy preservation

The nonlinear system ( 13) is solved with a modified Newton-Raphson method. The jacobian matrix is manually implemented and is only recomputed if the Newton iteration procedure does not lead to a satisfying solution. In order to decrease memory storage and the number of equations to solve, the degrees of freedom corresponding to the hammer and the shank as well as the corresponding equations are removed from the computation when the coupling force between the hammer head and the strings returns to zero for a long time.

In this paragraph, we only consider the hammer shank, head and piano strings. When the damping terms (included in the string, and in the hammer -string interaction) are set to zero, the global discrete energy of the { string, shank, head } system is conserved as shown in figure 4. This figure was obtained for note C♯5 (the hammer strikes 3 strings, see [START_REF] Chabassier | Physical parameters for piano modeling[END_REF] for more details on the parameters) which is the 53 rd key of a regular keyboard, using a constant jack force until let-off. The let-off distance was 3 mm for this simulation. When arithmetic with multiple precision is used for the numerical computations with a precision equal to 10 -60 , the energy seems to be very well conserved, confirming that the employed numerical method is conservative and stable. However, it can be observed some loss of accuracy, since the relative error on the energy is around 10 -54 . In double precision (the precision being equal to 2 • 10 -16 ), the relative error is approximately equal to 10 -8 . We think that an adimensionalization procedure should reduce this loss of accuracy.

Influence of pianistic touch

Our simulation tool enables us to study the influence of various parameters on the whole mechanical system. One of the most interesting parameters is the way the pianist strikes the key, also called "pianistic touch". As said earlier, the jack force is a data to our model, but it is only an intermediate force in the whole piano mechanism. Very few experimental recordings of this force are available and trustworthy (see [START_REF] Birkett | Experimental investigation of the piano hammer-string interaction[END_REF]), so we chose to use formula [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF] to represent the jack force. Moreover, we recall that the let-off mechanism is abruptly represented by the fact that this force is set to zero when the distance between the hammer head and the string becomes smaller than the let-off distance. The following results are obtained with note D♯1 which only has one string (see [START_REF] Chabassier | Physical parameters for piano modeling[END_REF] for more details on the parameters).

In the following, we will compare two simulations coming from the jack forces obtained with the parameters S and L listed in table 3. These two simulations have very different jack forces but the hammer head hits the string at about the same velocity (around 3.4 m/s), so that only the pianistic touch differs. They will schematically stand for staccato and legato touches, this classification will come from the comparison with the experimental measurements of figure 3.

Figure 5 shows different quantities obtained by simulation of note D♯1, using parameters S (left) and L (right) of table 3. In the upper frame S is displayed a global view of the string (length 1.965 m) and the hammer (shank and head). In the six frames (a) to (f) are displayed the shank and head deformations at different instants, in regular scale (black) and magnified scale (blue). In frames J and F are respectively displayed the jack force and the hammer -string interaction force with respect to time. The magnified shank deformation is clearly visible before the hammer -string 3 is very short and strong, while the jack force obtained with parameters L is progressively increasing. As a result, the jack force S is like a flick for the hammer shank, which is let free to move very soon (when the jack force returns to zero). On the contrary, the jack force L is still increasing when the shank reaches the let-off point, and therefore is suddenly set to zero because of this let-off mechanism. Comparing to figure 3 leads us to classify the jack force S as a staccato touch, and the jack force L as a legato touch.

This work was initially motivated by discrepancies between our experimental measurements and our simulated results. In the following, we will try to assess to what extent our new hammer shank -head model can diminish these dissimilarities, and evaluate the impact of the pianistic touch on the sound, if any.

Whereas it is very difficult to measure a jack force, it is rather easy to measure the shank acceleration. Putting a small accelerometer at the end of the shank enabled us to measure it on a Steinway D grand piano (see [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF] for more details and figure 6(b) for a photography of the experimental setup). For note D♯1 struck ff (very strong attack), we obtained the experimental signal plotted in figure 6(b). In article [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF], we noticed some "additional small oscillations" in the hammer acceleration which did not appear on the simulated hammer acceleration. At the time, no hammer shank was included in our piano model, the hammer was only represented by a hammer head given an initial velocity towards the strings, without accounting for gravity. We therefore made the hypothesis that these oscillations were due to the vibrations of the hammer shank. In figure 6(a), we have plotted the hammer head acceleration in the direction u θ , so that it can be directly compared to the experimental signal. The acceleration obtained with a legato touch is in blue and the acceleration obtained with a staccato touch is in red. As a reference, we have put on the same graph, in black, the hammer head acceleration obtained using the old model described in [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF], without any shank nor gravity, setting the initial hammer velocity close to the hammer head -string contact velocity of the other two simulations : v 0 = 3.41 m/s. A first remark is that both legato and staccato accelerations exhibit oscillations that the acceleration obtained with no shank does not. The acceleration simulated with the legato touch shows very fast oscillations which appear after let-off, because of the discontinuity created by this mechanism. These oscillations are very similar to the ones of the experimental signal, showing that accounting for the shank in the hammer -string interaction can explain the presence of strong oscillations in the hammer head acceleration. However, the simulated acceleration is clearly oscillating too much with the legato touch, most likely because our model of let-off is too abrupt compared to what happens in a real piano mechanism. A better model should account for the un-crushing of the roller felt (element of contact between the jack and the hammer shank), which is expected to make the jack force smoother during let-off.

A second remark is that different pianistic touches will lead to very different hammer head acceleration signals. An interesting question is to know whether this difference is noticeable as well in the hammer head -string interaction force, since this is the only transmitted quantity to the rest of the piano structure, and ultimately, the sound.

In figure 7 are displayed the hammer head -string interaction forces for legato and staccato touches (in blue and red) and with no shank (in black) with respect to time (figure 7(a)) and with respect to frequency (figure 7(b)). The three forces clearly differ in the time domain, but more spectacularly in the frequency domain, especially around 600 Hz where the staccato force exhibits a drop in its spectrum of 10 dB compared to the spectrum of the force obtained with no shank, and the legato force exhibits a drop even more dramatic of about 25 dB. The precise physical reason for this drop is still to understand and is probably linked to the nonlinear behavior of the shank. The reader should however be careful in reaching conclusions from these results because it is not physically relevant to compare the legato or staccato touches with the "no shank" case since the latter is purely idealistic and never found in practice. The only relevant comparison to do is between the different touches.

Simulations of the global piano model described in [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF] coupled with our new hammer model make it possible to compare the sound pressure signals obtained for these pianistic touches, and with no shank. Results are displayed in figure 8, where the first subfigure is obtained for a legato touch, the second one with a staccato touch and the last one with no shank. As expected, the observations made on the spectrum of the interaction forces are clearly visible in the sound pressure spectrum. The partials 14 to 19 are present at very different levels from one touch to the other. For instance, partial 16 is below -40 dB for the legato touch (let-off distance of 1.5 mm), but reaches almost -20 db in the staccato touch, and is above -20 dB when no shank is simulated. These discrepancies are audible in the isolated sounds for a trained ear (the sound coming from the staccato touch being brighter than the one coming from the legato touch), and the authors believe that they certainly induce a sound "color" when present in a musical phrase. Moreover, the touch seems to influence the levels of the longitudinal harmonics (these harmonics are responsible for the "bell" like color of bass piano sounds, and come from the propagation of waves in the elongation direction of the string): for this string, the fundamental frequency of the longitudinal wave is equal to 550 Hz. Its second harmonic, around 1100 Hz, and indicated by 2 f L in the figures, is clearly visible, audible, and at different levels from one pianistic touch to the other (-26.3 dB for the legato touch and -18.8 dB for the staccato touch, contributing to a brighter sound in the staccato case). To sum up, it is clear from these numerical experiments where all parameters, touch excepted, are rigorously the same, that the pianistic touch can remarkably influence the piano sound color. Notice that these spectral discrepancies only come from the differences in hammer shank vibrations. Indeed, they cannot come from other suggested causes in the literature as a shock excitation on the structure, nor from a changing impact point, nor a "longitudinal rubbing motion" on the string, since neither of these features are modeled in our study. The shock excitation on the structure measured in [START_REF] Koornhof | The influence of touch on piano sound[END_REF][START_REF] Goebl | Once again: The perception of piano touch and tone. Can touch audibly change piano sound independently of intensity ?[END_REF] is clearly associated to real staccato touches and its effects add up to the spectral change that we highlighted here thanks to our simulation tool.

Influence of the let-off distance

Another very interesting parameter to explore with our simulation tool is the let-off distance. This distance is a note by note setting done by the piano technician. We recall that in our model, as soon as the distance between the hammer head and the string (at rest) is smaller than the let-off distance, the jack force is set to zero. This is a crude modeling of reality but allows us to draw qualitative conclusions. Since the let-off distance is not relevant when using a staccato touch (because the force is already zero when the hammer head reaches the let-off point), the following tests have been made with three different let-off distances using a legato touch (simulations L, L+ and Lof table 3). Very extreme let-off distances have been tested : 3.0 mm, 1.5 mm as in the previous paragraph, and 0.0 mm, meaning that the jack force is set to zero at the exact contact time between the hammer head and the string. Interaction forces are displayed in figure 9 with respect to time (figure 9(a)) and with respect to frequency (figure 9(b)). There is a clear influence of the let-off distance on the interaction force, hence on the resulted sound. The difference on the force is rather spectacular in the [500 -1000] Hz region of the spectrum, with discrepancies up to 20 dB in very audible zones. The influence of the let-off distance on the resulted piano sound can be observed in figure 8, where the first subfigure is obtained for simulation L and the third one for simulation L+. The partials 14 to 19 are present at very different levels (see for instance partial 16 : 20 dB different), which has a great influence on the timbre.

Pianistic touch longitudinal precursor

In real piano mechanisms, the hammer head can have a planar motion (as opposed to a vertical motion only, in the model of article [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF]). Of course the main contribution of the hammer on the string is a vertical impulse, but there might exist a horizontal impulse coming from a rubbing motion of the hammer head on the string, as suggested in [START_REF] Askenfelt | From touch to string vibrations. II: The motion of the key and hammer[END_REF]. In our model, we have not considered this interaction, the system behaves as if the string was very slippery : the hammer head is free to have any horizontal motion, with no influence on the string. However, we can quantify this motion in order to estimate its extent on the resulted sound. To do so, two quantities are computed :

v θ = dξ dt • u θ and v r = dξ dt • u r
When the hammer and the string are in contact, the macroscopic angle θ is very small, thus v θ and v r can be interpreted as vertical and horizontal velocities. The vertical and horizontal velocities for legato and staccato touches are respectively displayed in figure 10(a) and 10(b) with respect to time. A first remark is that the horizontal velocity during contact (from 0 to 8 ms) is not small at all and oscillated between -1.5 and 1 m/s, which corroborates the hypothesis of a possible rubbing motion of the hammer head on the string. Another remark is that this horizontal velocity is quite different from one touch to the other, which leads us to believe that accounting for this rubbing motion could help us explain the influence of the pianistic touch in the piano sound. A first approach in order to model this rubbing motion would be to implement a continuity displacement relation in the horizontal direction between the hammer head and the string, as if they were not allowed to slip horizontally. The equations would have to respect an energy identity in order to be included in our model. Of course if such an interaction occurred, the hammer head would probably move less than what we can see in figure 10(b) where the horizontal motion is totally free. In such a model, the string would therefore be given an impulse in its elongation direction, making a "longitudinal" wave propagate in the string. The string model that we considered in [START_REF] Chabassier | Modeling and simulation of a grand piano[END_REF] already accounted for the presence of longitudinal waves (see figure 11(a)) because of their importance, via a nonlinear coupling with the transversal waves, in the generation of specific partials in the spectrum (so called "phantom partials") and in the presence of a nonlinear precursor in the time domain. This new model of rubbing would be a new way of generating longitudinal waves in the piano string that should significantly enrich the piano model and allow us to highlight the presence of a new precursor coming from the pianistic touch, and confirm the experimental intuition of [START_REF] Askenfelt | From touch to string vibrations. II: The motion of the key and hammer[END_REF]. Furthermore, this "pianistic touch longitudinal precursor" would certainly present a greater amplitude compared to the nonlinear precursor, which is only a second order phenomenon (see [START_REF] Chabassier | Energy preserving schemes for nonlinear Hamiltonian systems of wave equations: Application to the vibrating piano string[END_REF]) and would add up the one coming from the shock on the structure.

When the shank is modeled as a flexible beam, the hammer head -string contact point is not necessarily the top of the hammer head, because a rotation of the hammer head can occur, due to the bending of the hammer shank. Another improvement of the model related to the presence of a horizontal motion of the hammer head, would therefore be to compute the interaction force of equation ( 18) :

• with the real distance d between the hammer head and the string (no longer projected on the vertical axis),

• with coefficients p, K, and R that depend on the impact zone of the hammer head.

The first point is rather straightforward once noticed that the hammer head can rotate. The second point is less obvious and comes from the fact that piano technicians clearly divide the hammer head's felt into several zones that are provided with different cares during the "needling" operation. It consists in penetrating needles in the hammer felt (see figure 11(b)) in order to change its mechanical properties. In piano maintenance textbooks, three zones are identified, where the needling must be different (see figure 11(c)). Consequently, the coefficients of the hammer head -string interaction could be modeled as a function of the contact point on the hammer head's surface. Unfortunately, to the best of our knowledge, this kind of experimental data is not available in the literature and would be very difficult to calibrate manually.

Conclusions

In this paper, novel energy-preserving numerical schemes were provided for a Timoshenko beam in non-forced unknown rotation. Numerical simulations were carried out for the application case of a vibrating piano hammer shank coupled to a hammer head, coupled to a whole piano. To our knowledge, this work provides the first numerical simulation of such a complete system. Thanks to this simulation tool, we were able to highlight the effect of different pianistic touches, modeled with several input forces to the system (Jack forces). These numerical results show that the pianistic touch clearly influences the spectrum of the piano sound of equally loud isolated notes. These differences do not come from a possible shock excitation on the structure, nor from a changing impact point, nor a "longitudinal rubbing motion" on the string, since neither of these features are modeled in our study. Possible extensions of this work would be to refine the let-off model in order to make it smoother, and to consider the effect of hammer head's movements on the string: a possible transmission to the strings of the horizontal motion of the hammer head, a potential change of the impact point on the string, or even a change in the hammer stiffness coefficients if the contact point on the hammer's head surface were to vary.

Appendix B. Proof of consistency of system (12) with the inertial terms

We see that equation (12b) is not the natural consistent approximation of the continuous terms, but comes from the variational integration technique. Consistency is not obvious and can be shown as follows. We perform a Taylor expansion on the exact solution in order to evaluate the truncation error: Finally, we obtain the consistency of (12b):

θ n+1 -θ n ∆t = 1 ∆t ✟ ✟ ✟ θ(t n ) +
1 ∆t        θ n+1 -θ n ∆t w n+1 + w n 2 2 - θ n -θ n-1 ∆t w n + w n-1 2 2        = w 2 (t n ) θ(t n ) + 2 ẇ(t n ) w(t n ) θ(t n ) + O(∆t 2 )
Consistency of the discretisation (12a) is straightforward.

Figure 1 :

 1 Figure 1: Grand piano action. (a) Photography of an isolated Renner grand piano action, from kungfubrothers.com. (b) Sketch of the different parts of the grand piano action, from andersthorin.com. When the player hits the key, the jack pushes the hammer shank near its fixation point and forces it into a rotating motion.

  Schematic view of the mechanical system.

Figure 2 :

 2 Figure 2: Introduction of the notations.

FIG. 7 .FIG. 8 .FIG. 7 .Figure 3 :

 7873 FIG. 7. Registrations of acceleration at the hammer head and compression force in the jack (rnezzoforte, C4 ). The force was measured by means of apiezofoil in the jack. The position of the aeeelerometer is marked with a triangle. Also shown are the contact signals jack dolly and jack roller. (a) Legato touch and (b) staccato touch.

Figure 4 :

 4 Figure 4: Relative error on the energy versus time after let-off. Double precision (ε = 2 • 10 -16 ) and multiple precision (ε = 10 -60 ) are compared.

Table 3 :

 3 Parameters set t jack (s) A jack (N) let-off distance (m) Parameters of jack forces for different simulations contact (see frame (a) of figure 5(a)) but also during hammer -string contact (frames (c), (d) and (e) of figure 5(a) and 5(b)). The maximum time in frames J and F is different (25 ms for figure 5(a) and 45 ms for figure 5(b)) because the hammer does not reach the string at the same moment for both jack forces. The jack force obtained with parameters S of table

Figure 5 :

 5 Figure 5: Visualization of shank vibrations and hammer crushing after pianistic attack for staccato (left) and legato (right) touches. Frame S shows the whole string and the hammer at the same scale. The black rectangle is zoomed in at several times in frames (a) to (f). (a) : during jack force increase. (b) : at let-off distance. (c) : at hammer head -string contact. (d) : during the first decrease of the interaction force. (e) : at a second maximum of the interaction force. (f) : after hammer -string contact. Shank vibrations are magnified 30 times (in blue) and hammer crushing is magnified 7 times (in blue). Frames J and F respectively show the jack force and the interaction force with respect to time. The red vertical dashed lines correspond to the instants of frames (a) to (f).

  Experimental acceleration and photography of the accelerometer fixed on the shank

Figure 6 :

 6 Figure 6: Comparison between simulated (left) and experimental (right) acceleration of the hammer head. The time scale has been shifted so that time 0 corresponds to the contact instant between the hammer head and the string.

Figure 7 :

 7 Figure7: Temporal and spectral profiles of hammer head -string interaction forces obtained with different pianistic touches. The time scale has been shifted so that time 0 corresponds to the contact instant between the hammer head and the string.

Figure 8 :

 8 Figure 8: Spectral profiles between 0 and 1400 Hz of the sound pressure at point x = 0.850, y = 1.459, z = 0.380, in a cartesian coordinate system in which the soundboard is in the plane z = 0.200. Y-scales are in decibel (logarithmic) and normalized. Black dashed lines indicate the -20 dB and -40 dB levels in all the figures, for easy comparison.

  Power spectrum of the interaction forces

Figure 9 :

 9 Figure9: Interaction forces obtained with different let off distances (legato touch). The time scale has been shifted so that time 0 corresponds to the contact instant between the hammer head and the string.

Figure 10 :

 10 Figure 10: Velocity of the hammer head for Legato and Staccato touches. The time scale has been shifted so that time 0 corresponds to the contact instant between the hammer head and the string.

  Variables considered in the string model. Transversal displacement u s (m), longitudinal displacement v s (m), shear angle ϕ s (rad). (b) Hammer voicing with a needle voicing tool. From piesikpiano.com.

  Hammer head needling zones. Zone 1 : weak -light needling for piano dynamics. Zone 2 : medium -deep needling for mezzo-forte dynamics. Zone 3 : heavy -deep needling for forte dynamics. From[START_REF] Reblitz | Accord & Restauration, L'entretemps[END_REF].

Figure 11 :

 11 Figure 11: String model variables and theory of hammer voicing.

Table 1 :

 1 -) 0.086 10.18 • 10 9 32.38 • 10 -6 560 83.44 • 10 -12 0.64 • 10 9 0.85 Shank parameters

	H (m) d 0 (m) m H (kg) :	for D♯1	for C♯5
	0.04	0.02	10.76 • 10 -3 7.90 • 10 -3

Table 2 :

 2 Heads parameters

  2 (t n ) + ∆t ẇ(t n )w(t n ) + ∆t 2 4 ẇ2 (t n ) + 2w(t n ) ẅ(t n ) + O(∆t 3 ) = w 2 (t n ) θ(t n ) + ∆t 1 2 θ(t n )w 2 (t n ) + ẇ(t n ) w(t n ) θ(t n )

								∆t θ(t n ) +	∆t 2 2	θ(t n ) +	∆t 3 6	... θ (t n ) -✟ θ(t n ) + O(∆t 4 ) ✟ ✟
					= θ(t n ) +	∆t 2	θ(t n ) +	∆t 2 6	... θ (t n ) + O(∆t 3 )
		w n+1 + w n 2	2	= w(t n ) +	∆t 2	ẇ(t n ) +	∆t 2 4	ẅ(t n ) + O(∆t 3 )
	⇒	∆t θ n+1 -θ n	2 w n + w n-1		2
		+ ∆t 2 1 6	... θ (t n )w 2 (t n ) +	1 2	θ(t ∆t 2	θ(t n ) +	∆t 2 6	... θ (t n ) + O(∆t 3 )
		w n + w n-1 2 ẇ2 (t + ∆t 2 1 ∆t 2 4 6 ... θ (t n )w 2 (t n ) + 1 2 θ(t n ) ẇ(t n ) w(t n ) + 1 4 θ(t

2 = w n ) ẇ(t n ) w(t n ) + 1 4 θ(t n ) ẇ2 (t n ) + 2w(t n ) ẅ(t n ) + O(∆t 3 )

In the same way,

θ nθ n-1 ∆t = θ(t n ) -2 = w 2 (t n ) -∆t ẇ(t n )w(t n ) + n ) + 2w(t n ) ẅ(t n ) + O(∆t 3 ) ⇒ θ nθ n-1 ∆t w n + w n-1 2 2 = w 2 (t n ) θ(t n ) -∆t 1 2 θ(t n )w 2 (t n ) + ẇ(t n ) w(t n ) θ(t n ) n ) ẇ2 (t n ) + 2w(t n ) ẅ(t n ) + O(∆t 3 )

for Timoshenko systems, it is known that locking phenomena can appear (see[START_REF] Chapelle | An optimal low-order locking-free finite element method for Reissner-Mindlin plates[END_REF]), but using high order Finite Element Method diminishes this effect (see[START_REF] Cohen | Higher Order Spectral Finite Elements for Reissner-Mindlin Equations in the Time Domain[END_REF]).

Acknowledgments

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

Appendix A. Derivation of the model : virtual work principle

In this appendix we explain all terms involved in the virtual work identity.

Appendix A.1. Acceleration

The virtual work coming from acceleration has two contributions : from the shank and from the head. The calculus is completed from the expression of the second-derivative of r: r = -τ φ cos ϕ + τ φ2 sin ϕ -(s -τ sin ϕ) θ2 + 2( ẇ -τ φ sin ϕ) θ + (w + τ cos ϕ) θ u r + 2 τ φ cos ϕ θ -(s -τ sin ϕ) θ + ẅ -τ φ sin ϕτ φ2 cos ϕ -(w + τ cos ϕ) θ2 u θ and the virtual displacement δr:

The values of ϕ are assumed to be very small (Timoshenko hypotheses), the following substitutions are performed :

Therefore, the second-derivative of r is approximated as:

The integral of 1 on the section of the beam is equal to A the area of this section, the integral of τ equal to 0 (due to symmetry) and the integral of τ 2 is equal to I the bending moment of inertia of the beam. For instance, in the case of a circle cross-section, we have

, where R is the radius.

20

The virtual work due to the beam acceleration is then equal to:

and the virtual work due to the head acceleration is

. Exterior forces on the shank

The gravity field applies at each point of the beam:

So the virtual work of gravity on the beam reads:

The jack force is located at point x = x jack and oriented along -u θ :

The coupling force with the head is located on the center of gravity of the head and acts on the last point of the shank. It has to be written in terms of shank's unknowns:

. Exterior forces on the head

Gravity applies to the head's mass:

Furthermore, the opposite of the coupling force acts back on the head:

We take into account the force applied by the string:

F string is considered as a given force amplitude.

Appendix A.4. Interior forces

The interior forces of the beam come from the elastic deformation. Under the hypothesis of a linear Hooke's law, and following the literature (see [START_REF] Lin | Vibration analysis of a rotating Timoshenko beam[END_REF]), we write: