A Ozel

J C Brändle De Motta

Micheline Abbas

Pascal Fede

Olivier C Masbernat

Stéphane Vincent

J.-L C Estivalezes

O C Simonin

Particle resolved direct numerical simulation of a liquid-solid fluidized bed: Comparison with experimental data

Keywords: Particle-resolved DNS, liquid-solid fluidized bed, particle and fluid agitation

Introduction

Liquid fluidization is used in various industrial application involving biochemical, catalytic reactions and crystallization processes. The flow in a liquid fluidized bed lies within an intermediate regime between the settling of particles controlled by the hydrodynamic interactions and the rapid granular flow controlled by the collisions between particles, where the particle Reynolds number is in a range of O(100) and the particle Stokes number is in a range of O(10), both based on particle settling velocity. In this sense, liquid fluidization is a challenging problem for two-phase modeling. For practical applications, twophase continuum models are generally used to carry out numerical simulations, based upon two-fluid or statistical models [START_REF] Gevrin | Granular pressure and particle velocity fluctuations prediction in liquid fluidized beds[END_REF][START_REF] Zhang | Two-and three-dimensional computaional studies of liquid-solid fluidization[END_REF]. However, modeling of liquid-solid fluidization is still an open research topic and multi-scale modeling developments are still needed to correctly predict inter-particle and particle-fluid interactions. One major issue is to predict the right level of particulate and carrying flow phase fluctuations as a function of bed solid phase fraction (or fluidization velocity).

Resolved particle direct numerical simulations of particulate flows has been developing last two decades (see the review of [START_REF] Tenetti | Particle-resolved direct numerical simulation for gas-solid flow model development[END_REF]).

These simulations can provide the particulate phase fluctuation characteristics in order to develop appropriate two-phase continuum models. Many of particle resolved simulations have been carried out on fixed structured grids to take advantage of parallelization and avoid the complexity of mesh reconstruction. [START_REF] Pan | Fluidization of 1204 spheres: simulation and experiment[END_REF] carried out resolved simulations of fluidization of 1204 finite size spheres in a 2-D bed using the method of distributed Lagrange multipliers and as simulation results, the fluidization velocity versus fluid fraction was found to be a power law which exponent well compared with that predicted by the correlation of [START_REF] Richardson | The sedimentation of a suspension of uniform spheres under conditions of viscous flow[END_REF]. [START_REF] Zhang | Microstructural effects in a fullyresolved simulation of 1,024 sedimenting spheres[END_REF] performed a 3-D fully resolved simulation of 1024 particles settling under gravity in a periodic domain accounting for elastic collisions of particles. Their method is based on a linearization of Navier-Stokes equations in the vicinity of particle interface [START_REF] Zhang | A second-order method for three-dimensional particle simulation[END_REF]). In their study, Particle Reynolds number and solid volume fraction were respectively 10% and 13%. They have shown that the settling velocity was matching [START_REF] Richardson | The sedimentation of a suspension of uniform spheres under conditions of viscous flow[END_REF] correlation and evidenced the relation between the velocity fluctuations and particles microstructuration.

Using a Lattice Boltzmann Method to solve the interstitial flow and an equation of motion accounting for lubrication and collisions between particles, [START_REF] Derksen | Direct numerical simulations of dense suspensions: wave instabilities in liquid-fluidized beds[END_REF] have simulated in limited size domains the propagation of concentration waves in liquid-solid fluidized beds with large bed solid fraction (close to maximum packing) and particle Reynolds number of order of O(10). Their results were in qualitative agreement with an experimental study of [START_REF] Duru | Experimental investigation on the secondary instability of liquid-fluidized beds and the formation of bubbles[END_REF]. Based on the same method, [START_REF] Derksen | Simulations of scalar dispersion in fluidized solid-liquid suspensions[END_REF] performed the simulation of the mixing of a passive scalar in a fluidized bed with periodical boundaries in a wide range of bed solid volume fraction (0.2-0.5) and particle Reynolds numbers of order 10. Derksen's results first show a good agreement with [START_REF] Richardson | The sedimentation of a suspension of uniform spheres under conditions of viscous flow[END_REF] exponent dependence with Reynolds number. Interestingly, [START_REF] Derksen | Simulations of scalar dispersion in fluidized solid-liquid suspensions[END_REF] showed that the diffusion of the passive scalar in the bed is similar to the auto-diffusion of particles, scaling of which is close to what was experimentally observed in sedimentation by [START_REF] Nicolai | Particle velocity fluctuations and hydrodynamic self-diffusion of sedimenting non-brownian spheres[END_REF]. [START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] developed an Immersed Boundary Method to simulate the sedimentation of 1000 spherical finite particles at high Reynolds number (400) and highly dilute limit, but no quantitative comparison with existing data was provided. More recently, [START_REF] Uhlmann | The motion of a single heavy sphere in ambient fluid: A benchmark for interface-resolved particulate flow simulations with significant relative velocities[END_REF] evaluated the accuracy of their method as a function of the spatial resolution (number of meshes per particle diameter) for the case of a single sphere settling in an infinite stagnant fluid, in a wide range of Reynolds and Archimedes (or Galileo) numbers. The higher the latter number, the higher spatial resolution is required, up to 48 mesh points per particle diameter at high Galileo number. Then [START_REF] Chouippe | Forcing homogeneous turbulence in direct numerical simulation of particulate flow with interface resolution and gravity[END_REF] used this method to study turbulent particle settling in a channel. [START_REF] Corre | Direct numerical simulation of a liquid-solid fluidized bed computational techniques for multiphase flows[END_REF] used a fictitious domain approach to perform particleresolved simulations of the liquid-fluidized bed experimentally studied by [START_REF] Aguilar-Corona | Agitation of particles in a liquid fluidize bed[END_REF]. Instantaneous and averaged flow characteristics of the fluidized bed were qualitatively in good agreement with experimental trends. Since then, this method was improved and has been applied in the present study with a higher level of accuracy [START_REF] Vincent | A Lagrangian VOF tensorial penalty method for the DNS of resolved particle-laden flows[END_REF]). The numerical technique is a four-way coupling method, based on a one-fluid formulation of the incompressible Navier-Stokes equations solved on a structured Cartesian grid.

Flow parameters

Flow parameters chosen for the simulation of the fluidized bed are taken from the experimental study of [START_REF] Aguilar-Corona | Agitation of particles in a liquid fluidize bed[END_REF] the drag coefficient of a single particle at V t , here equal to 0.8). Fluidization law and fluctuating motion of both phases have been measured by [START_REF] Aguilar-Corona | Agitation of particles in a liquid fluidize bed[END_REF] in a range of fluidization velocities ranging between 0.17 and 0.05 m/s, corresponding respectively to be solid volume fraction ranging between 0.1 and 0.5. Details of the measurement techniques can be found in [START_REF] Aguilar-Corona | Agitation of particles in a liquid fluidize bed[END_REF] et al. (2014). The DNS approach is based on a one-fluid formalism of the incompressible Navier-Stokes equations with an algebraic adaptive augmented Lagrangian method used for pressure-velocity coupling (Implicit Tensorial Penalty Method, ITPM). The particles are considered as a fluid with specific rheological properties whose evolutions are modeled by the Navier-Stokes equations. This method enforces the solid behavior of the particles in the framework of Eulerian fixed grid. A Lagrangian Volume Of Fluid (VoF-Lag) method enables particle tracking while avoiding particle shape deformation and ensuring volume conservation of the solid phase. This approach provides a second order convergence in space and time.

U F 0 (1 -φ b ) n n = 2.41, U F 0 = 0.226 m/s
In dense flows such as in fluidized beds, wall-particle and inter-particle collisions must be accounted for as well as lubrication. Fully resolved lubrication fluid layer between a particle and a wall requires a highly refined Eulerian grid at the scale of the particle (about 150 grid points per particle that well matches the scaling law proposed by [START_REF] Legendre | A note on the modelling of the bouncing of spherical drops or solid spheres on a wall in viscous fluid[END_REF]. Using the parameters of the studied fluidized bed and varying the normal collision velocity, the resulting normal coefficient of restitution is also in good agreement with the correlation of [START_REF] Legendre | A note on the modelling of the bouncing of spherical drops or solid spheres on a wall in viscous fluid[END_REF], as illustrated by Figure 1. The collision and lubrication force models for multiple particle-particle and particle-wall interactions are implemented in the Navier-Stokes equations as volume force terms, and their semi-implicit treatment avoids particle overlapping during the solving step of the flow field.

Figure 2 shows the computational domain. It has the same dimensions as the experimental fluidization column. It is composed of a parallelepipedic box of dimensions 0.08 × 0.08 × 0.64 m. The solid wall boundary of the cylindrical column is simulated using a Darcy penalty method [START_REF] Khadra | Fictitious domain approach for numerical modelling of Navier-Stokes equations[END_REF], consisting in adding a Darcy term in the momentum equations with a very small permeability ascribed to the cells located outside the cylindrical envelope. This method ensures a no-slip condition at the cylinder wall. A uniform distribution of fluid velocity is imposed at the bottom of the bed and a free outlet boundary condition is defined at the top of the bed. The domain is discretized with Statistical averages presented in next section were therefore calculated over a 10 s period of simulation, 10 s after the beginning of each simulation (see figure 4).

Results and Discussions

Numerical results are compared to experimental data taken from [START_REF] Aguilar-Corona | Agitation of particles in a liquid fluidize bed[END_REF]. In experiments, averaged bed height was determined from video camera in slightly unmatched refractive index conditions. Measurement of this quantity at different fluidization velocities gives the fluidization law. Fluctuating motion of the particles was analyzed from the recording of 12 trajectories of marked particles during 3 minutes at a sampling frequency of 60 Hz. Fluctuating motion of the liquid phase was characterized from the acquisition of the velocity field in a median plane of the column, using high speed PIV (between 250 and 500 Hz) and a spatial resolution of d p /5.

Fluidization law

The fluidization law (relationship between the fluidization velocity and bed solid fraction) is the first step of validation of the numerical model, reflecting the macroscopic balance between buoyancy and drag forces. In order to calculate the bed solid concentration φ b , the bed height h b was computed using two different methods, which led to the same result.

First, the bed height was set equal to the time averaged maximum particle position in the axial direction. U F = 0.17, 0.15, 0.12, 0.09, 0.073 m/s (resp. φ b = 0.11, 0.14, 0.22, 0.31, 0.39)

U F = 0.17 [m/s] U F = 0.15 [m/s] U F = 0.12 [m/s] U F = 0.09 [m/s] U F = 0.073 [m/s]

Getting statistics

φ b = 2 3 n p d 3 p h b D 2 (1)
Plotting the fluidization velocity as a function of the bed solid concentration φ b gives the fluidization law. Figure 6 shows the bed solid concentration measured in DNS together with the experimental data of [START_REF] Aguilar-Corona | Agitation of particles in a liquid fluidize bed[END_REF] and the correlation of [START_REF] Richardson | The sedimentation of a suspension of uniform spheres under conditions of viscous flow[END_REF]:

U F = U F 0 (1 -φ b ) n (2)
where n is a function of Re t , U F is the fluidization velocity and U F 0 is the fluidization velocity leading to particles entrainment. The exponent value best fitting experimental data is n = 2.41 (Re t = 530), and is in quite good agreement with the value predicted by Richardson-Zaki correlation at that particle Reynolds number (n = 2.39 for Re t > 500). Experimental value of U F 0 is found equal to 0.226 m/s, and the measured terminal velocity of the particles is also that it is not related to the ratio d p /D (Di Felice and Kehlenbeck ( 2000)).

V t = 0.24 m/s. The ratio of U F 0 /V t is equal to 0.
The agreement between experiments and numerical simulations is quite good, and to the best of our knowledge, this result is the first validation of particle resolved simulation of fluidization law in a full 3-D fluidized bed in that range of particle Reynolds number. becomes more and more angular as the bed is compacted, in response to the increase of inter-particle collisions. Overall, the multi-scale diffusive-like motion of the particles in the bed as calculated by the numerical simulations exhibits remarkable similarities with the experimental signals (12 particle trajectories recorded during more than 3 minutes), suggesting that the physics of the fluidparticle and particle-particle interactions are qualitatively well captured by the numerical model, in both dilute and dense regimes. the collisional model (used in the numerical simulations) that generates abrupt modification of the velocity upon collisions, whereas their absence in the experimental signal could be due to a filtering effect (velocity sampling frequency in the experiments: 30 Hz).

Particle trajectories and Lagrangian velocity signal

Recirculation

The low frequency, large amplitude fluctuations suggest the presence of largescale coherent structures. cylinder is weak, and statistics of particle velocity are not fully converged near the bed axis in the 10 seconds integration window. This lack of convergence near the bed axis is more pronounced at highest fluidization velocities (lowest bed solid fraction). As a general trend particle motion in the bed is upward in the middle of the bed and downward in the near wall region (between 0.2 and 0.6 column radius from the wall). When particles are very close to the wall they tend to rise along the wall. Note that the magnitude of this mean motion is an order of magnitude smaller than the particle r.m.s. velocity presented in the next section.

Particle and fluid velocity variance

The average of the velocity variance of particles in the whole bed is computed as follows:

< u 2 p,i >=< (u p,i -< u p,i > ins ) 2 > (3)
< . > denotes the time-space average defined in the Appendix, u p,i is the i component of the instantaneous Lagrangian particle velocity, and the symbol and z/h b = 0.9. In the bottom part of the bed, 0 < z/h b < 0.2, the variance is smaller and is growing from 0.5 to 0.9. In the top part of the bed (z/h b > 0.9), the occurrence of peaks of large amplitude in the near-freeboard region results from the transition of fluid velocity between a concentrated medium and a free- of the profile is rather homogeneous varying from 0.9 in the core of the bed to 1.1 near the wall.

Figure 13 shows the variance of particle velocity components in the whole bed (eq. ( 3)) as a function of bed solid volume fraction and compared to experimental data. The numerical data are reasonably in good agreement with experiments, showing a strong decrease of the agitation with the concentration.

However at the largest fluidization velocity (lowest concentration), the simulation overpredicts the experimental value by a factor of two. At larger concentrations the trend is reversed, numerical results underpredict experimental values.

The same behavior is observed for both components (axial and transverse).

The particle agitation in the fluidized bed is not isotropic: the variance of the axial component of velocity is stronger than that of the transverse component (isotropy of the fluctuations in the transverse plane was checked). This behavior is already well known in gravity driven gas-solid suspensions, like in sedimentation. The anisotropy coefficient k anis is defined as the ratio of the particle velocity variance in flow direction (z) to that in the transverse plane (x, y):

k anis = < u 2 p,z > 1 2 (< u 2 p,x > + < u 2 p,y >) . ( 4
)
Its evolution as a function of bed solid fraction is shown in Figure 14. Numerical predictions (1.5 in average) are close to the experimental values (1.6) and nearly constant in the range of bed solid fraction investigated, unlike sedimentation case at low Reynolds number where anisotropy is decreasing as particle concentration increases [START_REF] Nicolai | Particle velocity fluctuations and hydrodynamic self-diffusion of sedimenting non-brownian spheres[END_REF]). This result emphasizes the leading role of large-scale motion in the agitation of particles fluidized by a liquid at high Re p . Note that despite the differences observed between numerical and experimental data at lowest and highest solid fraction, the correct prediction of anisotropy coefficient over all the concentration range investigated suggests that the structure of the large scale motion is well captured by the numerical model.

The total fluid energy, E f and particle fluctuating kinetic energy,E p with respect to the bed solid volume fraction are shown in Figure 15. For the particle phase it is computed as:

0 0.2 0.4 0.6 0.8 1 0 0.5 1 1.5 2 U F = 0.17 m/s U F = 0.15 m/s U F = 0.12 m/s U F = 0.09 m/s U F = 0.073 m/s < u 2 p,x + u 2 p,y > layer / < u 2 p,x + u 2 p,
E p = 1 2 2 u 2 p,x + u 2 p,z (5)
and for the fluid phase:

E f = 1 2 (2 u 2 f,x plane + u 2 f,z plane ) (6)
where [.] plane denotes the Eulerian average in a vertical median plane (see definition in Appendix). The choice of this average is driven by the correspondence with experimental data obtained with High Frequency PIV in a vertical median plane of the bed (Aguilar-Corona ( 2008)). The particle and fluid agitation are shown here at moderate concentration range, while one can expect two limiting behaviors at small and high volume fractions.

φ b E p , E f [m 2 /s 2 ] Sim., E p Sim., E f Exp., E p Exp., E f
In the very dilute regime (φ b → 0), particle agitation is expected to be close to that observed in a dilute turbulent pipe flow. This regime corresponds to a fluidization velocity equal to the particle terminal velocity (particle entrainment). In the present case, particle terminal velocity is V t = 0.24 m/s, and flow

Reynolds number is Re

f = ρ f VtD µ f ∼ 7×10 3
, with a wall friction velocity derived from Blasius law equal to 0.016 m/s. The corresponding fluctuating kinetic energy is about 3 × 10 -4 m 2 /s 2 , which is nearly two orders of magnitude smaller than E f measured (in the experiments and simulations) at φ b = 0.1. Neglecting the fluid turbulent modulation by the particles in the very dilute limit, an estimation of particle fluctuating kinetic energy can be scaled as that of the continuous phase weighted by a function of particle Stokes number [START_REF] Tchen | Mean value and correlation problems connected with the motion of small particles suspended in a turbulent fluid[END_REF], [START_REF] Deutsch | Large eddy simulation applied to the motion of particles in stationary homogeneous fluid turbulence[END_REF]), here defined as the ratio of particle response time to the fluid turbulent time macro-scale. For the present system, such an estimate gives in very dilute regime the same order of magnitude for and the other one due to small scale fluid-particle wake turbulence (also referred to in the literature as pseudo-turbulent kinetic energy), E f = E f + δE f , the second component would be far greater than the first one. Note that a rigorous formulation of the energy decomposition is reviewed in [START_REF] Fox | On multiphase turbulence models for collisional fluidparticle flows[END_REF]. The strong coupling seems to be intrinsic to liquid fluidization whereas the opposite is true in gas fluidization, which is essentially related to difference in particle inertia.

E p | φ b →0 as E f | φ b →0 (taken equal
When the concentration increases, the particle fluctuating energy decreases.

At large volume fraction (φ b → φ max ), the particle phase approaches a porous media. E p vanishes whereas E f remains finite, meaning that velocity fluctuations of both phases become uncorrelated. In this limit, large scale motion disappears and flow fluctuations derive from the so-called pseudo-turbulence (E f = δE f , or equivalently E f ≈ 0). Interestingly, the decrease of fluid fluctuating energy with concentration

E f | φ b =0.1 -E f | φ b →φmax is close to E p | φ b =0.1 ,
suggesting that the particle fluctuating energy at low concentration is mainly driven by the flow large scale fluctuating motion E f .

Note that in figure 15, the decay with bed solid fraction of particle agitation is stiffer in numerical than in experimental curves, whereas the reverse trend is observed with fluid agitation. The origin of the differences observed between numerical and experimental data is difficult to identify. First, statistics on the particle phase are not derived in the same way (2133 particles during 10 seconds for the numerical data, 12 trajectories during 3 minutes for the experimental data). Second, in this range of particle Reynolds number (Re t = 530), the flow is probably under-resolved with 12 cells per particle diameter and the smallscale structures of the flow are probably partially filtered [START_REF] Uhlmann | The motion of a single heavy sphere in ambient fluid: A benchmark for interface-resolved particulate flow simulations with significant relative velocities[END_REF]). The resulting particle relative velocity prior to collisions, and therefore the numerical treatment of collisions can be affected. All these issues require to be addressed separately in order to quantify their contribution to the calculation of both phases agitation.

Particle fluctuation time scales

The time-scale (macro-scale) characteristic of particle agitation can be derived from the computation of autocorrelation function of axial and radial velocity components (given in Appendix). The autocorrelation function is shown in Figure 16 for two bed solid volume fractions (0.12 and 0.3, respectively corresponding to 0.17 and 0.09 m/s). After the initial step of continuous decay, both components exhibit a large oscillating behavior at long times. Origin of these oscillations is likely due to the contribution of the localized recirculation zones in the bed bottom section (see Figure 10), which seems to be supported by the observation that period of oscillations decreases as the fluidization velocity is decreasing, i.e. when the bed height is decreasing. Comparison of numerical and experimental curve shapes exhibits some discrepancies. At short times (< 0.4 s) curves behave the same with a stronger decay of the transverse component compared to the axial velocity component. At longer times, experimental curves decay much more slowly than calculated curves, which oscillate around the abscissa axis. This is particularly noticeable for the autocorrelation function of the particle axial velocity component. Regarding the radial motion, the decorrelation time of velocity fluctuations, that is to say the time at which the curves reach the horizontal axis is nearly constant for both bed solid fractions and compares well with experimental data, close to 0.3 s. The decorrelation time of axial fluctuations as predicted by the numerical data is significantly smaller compared to the experimental data and its evolution with the bed solid fraction is reversed. When the bed solid volume fraction increases from 0.12 to 0.3, it decreases from 0.5 s to 0.3 s whereas experimental data shows an increase from . restitution coefficient becomes smaller than 0.05 (see Figure 1), which makes the collisions difficult to detect from either numerical or experimental signals.

Below 0.3, the dimensionless collision frequency is well predicted by the numerical model. It remarkably fits the theoretical expression derived from the kinetic theory (noted KTFG in Figure 17). From the scaling of the collision time,

d p / 3 4 u 2 p,x + u 2 p,y
that can be identified to d p / 3 2 θ p in KTGF, one can conclude that the transverse fluctuating motion of particles is the correct characteristic velocity scale to be considered for the collisions in the fluidized bed. This is consistent with [START_REF] Février | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF] and Fox (2014) who suggested that the total particle velocity fluctuations can be decomposed in large and small scale fluctuations E p = E p + δE p . The first part contains particle large scale motion represented by the streamlines of Figure 10. It is approximately equal to the particle velocity variance in the axial direction and is fully coupled to the flow large scale motion via the buoyancy force and the non-uniform two-phase mixture density field ( E p ≈ E f ). The second part accounts for random uncorrelated motion (transverse fluctuations), similar to Brownian motion resulting from collisions, usually referred to as granular temperature 3 2 θ p in gas-solid flows. Note that the measured and calculated particle velocity variances (figure 13) suggest that E p /δE p > 1 in liquid fluidization but this ratio is expected to decrease with particle inertia [START_REF] Février | Partitioning of particle velocities in gassolid turbulent flows into a continuous field and a spatially uncorrelated random distribution: theoretical formalism and numerical study[END_REF].

Conclusions

Particle resolved simulations of a liquid-solid fluidized bed were performed using a one-fluid formulation of the incompressible Navier-Stokes equations, where the pressure-velocity coupling is provided by an algebraic augmented Lagrangian method and particles presence is modeled with an implicit penalty fictitious domain method, sub-grid scale lubrication force and soft-sphere collision models. We carried out simulations in a fully 3-D fluidized bed experimentally investigated by [START_REF] Aguilar-Corona | Agitation of particles in a liquid fluidize bed[END_REF] on a structured uniform Eulerian grid at

  The resolved-scale particles are modeled by an Implicit Tensorial Penalty Fictitious Domain Method (ITPM). They are tracked by using a hybrid Eulerian-Lagrangian Volume of Fluid approach, which accounts for collisions and lubrication effects. This study has two scopes. The first one is to evaluate the effective ability of ITPM to predict two-phase flow behaviour by performing particle resolved simulations of a liquid-solid fluidized bed involving finite size particles, with large particle Reynolds and moderate Stokes numbers. The second one is to analyse velocity fluctuations of both phases in this regime. The bed geometry, particle size and number and flow parameters used in these simulations are the same as in Aguilar's experiments, allowing a direct quantitative comparison between experiments and numerical data. The paper is structured as follows: Flow parameters and numerical model (detailed in other references) are briefly presented in sections 2 and 3 respectively. Statistical quantities (as defined in appendix A) computed from the numerical results are compared with experimental data obtained by Aguilar Corona with same flow parameters and geometry. Fluidization law and particle velocity fluctuations predicted by the simulations are also compared in section 4.

  in a cylindrical column of 8 cm inner diameter. Phase material properties and fluidization parameters are reported in Table 1. Monodisperse spherical beads of Pyrex (d p = 6 mm, ρ p = 2230 kg/m 3 ) have been fluidized in a concentrated aqueous solution (65% w/w) of potassium thiocyanate (KSCN) of density ρ f = 1400 kg/m 3 and viscosity µ f = 3.8 × 10 -3 P a.s at T = 20 • C. At this temperature, refractive indices of both phases are matched, allowing the implementation of optical techniques such as high-speed video for the 3-D Lagrangian tracking of colored particles or high frequency Particle Image Velocimetry for the measurement of the velocity field in the liquid phase (Aguilar-Corona, 2008). Particle terminal velocity, V t , is 0.24 m/s and Reynolds number based on V t is Re t = 530. Inertia of the particles is characterized by a Stokes number here defined as St t = 8ρp 3ρ f C Dt = 5.3 (C Dt is

Figure 1 :9

 1 Figure 1: Normal restitution coefficient of collisions with respect to collisional Stokes number St coll = 2 9 Ru coll ρp µ f : influence of the spatial resolution.

Figure 2 :

 2 Figure 2: Geometry and dimensions of the computational domain.

Figure 3 :

 3 Figure 3: 3-D views of the simulated particles inside the bed (from left to right: t = 0, 5, 10 and 20 s).

  Figure 4 displays the time-evolution of this parameter at different fluidization velocities. After a transient period, it oscillates around a steady value for all cases investigated due to particle agitation. The intensity of fluctuations is a decreasing function of the fluidization velocity, with a maximum of the order of 5% for the lowest fluidization velocity. Second, the bed height was also determined by averaging in time (during 10 s) and space (over the bed volume) a particle phase indicator function χ p (x, t). This function is defined on each Eulerian mesh cell, equal to 1 if the node is inside the particle and 0 if not. A vertical profile of the time-section average of the particle phase indicator function (or phase fraction) {χ p } layer as defined in equation (A.3), is shown by Figure 5 at different fluidization velocities. This quantity represents the solid volume fraction averaged in cylinders of diameter D and thickness ∆z. So the integral of {χ p } layer layer along z is equal to φ b h b .

Figure 5

 5 Figure5shows that the phase fraction is rather homogeneous along the bed height but in the freeboard region a gradient of the solid volume fraction develops, becoming stiffer as the fluidization velocity is decreased. The bed height is then computed by applying a linear regression with a high order polynomial interpolation on the volume fraction profile at the interface between the bed and the freeboard region, and the z value of the inflection point of that function is defined as the bed height h b .Both estimations lead to close values of h b for all velocities with a difference of a few percent. The bed solid fraction φ b is then calculated according to:

Figure 4 :

 4 Figure 4: Time-evolution of maximum axial particle position at different fluidization velocities:

  Figure 5: Axial profile of time-section average of particle phase fraction χp.

Figures 7 Figure 6 :

 76 Figures 7 and 8 exhibit projections of 16 particle trajectories in the radial and vertical planes of the bed, for two different fluidization velocities (0.15 and 0.073 m/s). It can be observed that particle trajectories are quite sensitive to this parameter. At high fluidization velocity (low concentration), trajectories occupy all the bed space, with an apparent slight deficit of particles in the bed bottom zone (close to the flow inlet). For the same simulation time (10 s), the space travelled by the same number of particle trajectories tends to reduce at higher concentration (lower fluidization velocity). This confinement effect can be clearly observed on the trajectory envelopes projected in the cross section (x-y plane), with the development of dark spots near the bed wall, the signature of particle trapping over long-time periods. Additionally, the shape of the paths

Figure 9 Figure 7 :Figure 8 :Figure 9 :

 9789 Figure9exhibits the instantaneous axial and radial particle velocity components following one particle trajectory for the case φ b = 0.31 (U F = 0.09 m/s).Numerical and experimental signals present qualitatively similar features, being composed of large scale, low frequency and small scale, higher frequency fluctuations. The amplitude of fluctuations is more pronounced on the axial component U p,z than on the transverse one, U p,x . The frequency of the high amplitude velocity fluctuations is smaller on U p,z than on U p,x signals.Modes of high frequency can be observed in the numerical signal, which is not the case of the experimental signal. Their occurrence could be inferred to

  Figure10shows time and azimuthal average of the particle velocity field defined as u p ann∆z . It represents the average of particle velocity in a hollow cylinder of inner radius r and outer radius r + ∆r and of thickness ∆z (c.f. Appendix). Fields of u p ann reveal the presence of a localized large-scale recirculation in the lower section of the bed, size of which compares with the bed diameter. For both concentrations φ b = 0.31 (U F = 0.09 m/s) and φ b = 0.39 (U F = 0.073 m/s), particles preferentially rise up near the bed axis and flow downward near the wall. The upward velocity is larger than the downward velocity owing to mass conservation. The shape of particle path-lines indicates that the recirculation is stronger in the bottom of the bed, confirming the existence of a large-scale toroidal motion above the flow inlet detected in the trajectography experiments (cf Figure7). Figure11displays radial profiles of the axial component of particle velocity averaged in time and over the bed height u p,z ann ( . ann denotes the average of . ann∆z over bed height). When radial position tends towards zero, the density of particles in the control hollow

<

  . > ins refers to the average operator as < . > at each time step. The variance of radial particle velocity as a function of vertical position in the bed, calculated in horizontal layers of thickness ∆z and diameter D,is noted u 2 p,x + u 2 p,y layer(see Appendix). The axial profile of this quantity normalized by the radial velocity variance in the whole bed is reported in Figure12(top) at different fluidization velocities. Profiles collapse on a single curve, and exhibit a rather homogeneous distribution of radial velocity variance in a large portion of the bed, slightly increasing from 0.9 to 1.1 between z/h b = 0

Figure 10 :Figure 11 :

 1011 Figure 10: Averaged particle velocity field and streamlines for the case U F = 0.09 m/s (top) and U F = 0.073 m/s (bottom).

  particle domain, getting sharper as the bed solid fraction increases (or as the fluidization velocity decreases). The variance of axial particle velocity component as a function of radial position in the bed is calculated in vertical hollow cylinder of thickness ∆r and is noted u 2 p,z ann (see Appendix). The radial profiles of normalized variance figure 12 (bottom) for various fluidization velocities. The axial velocity variance is minimum in the core of the bed and maximum near the wall, due to the transition between negative velocities in the recirculation loop to positive values very near the wall (see Figure 11). However, the shape of the profiles depends on the fluidization velocity. At highest fluidization velocity (minimum bed solid fraction), the profile shows a marked gradient along bed radius (from 0.55 to 1.25), which tends to flatten as the fluidization velocity decreases. At largest bed concentration (smallest fluidization velocity) the shape

Figure 13 :Figure 14 :

 1314 Figure 12: Top: Axial profile of transverse particle agitation; Bottom: Radial profile of axial particle agitation.

Figure 15

 15 Figure15shows that the fluid fluctuating kinetic energy with respect to the solid volume fraction does also fit well experimental data obtained from HF PIV measurements, except at low phase fraction where it underestimates experimental data. Interestingly, the fluctuation level of the liquid phase is always significantly larger than that of the solid phase in all the range of solid phase fraction investigated. The ratio E f /E p is a growing function of φ b .

Figure 15 :

 15 Figure 15: Fluctuating kinetic energy of fluid phase E f and particle phase Ep. Solid symbols: simulations, empty symbols: Aguilar-Corona (2008).

  to that of a steady turbulent pipe flow), and therefore is one order of magnitude smaller than the particle fluctuating kinetic energy measured at the lowest volume fraction E p | φ b =0.1 . Therefore, both particle and fluid agitation should increase at very low concentration, without however being captured by the present measurements. From these estimations, we conclude that the fluctuating energy for both phases at φ b = 0.1 is far larger than the very dilute limit, which means that the fluctuations are dominated by strong particle-fluid flow coupling. Consequently, if we decompose the fluid fluctuating energy in two components, one of them induced by large scale collective motion

1. 8 Figure 16 : 24 Figure 16 :

 8162416 Figure 16: Autocorrelation fonction. Top : experiments. Circles, axial component, tirangles, radial component. Black : b = 0.3, Red= b = 0.12 Bottom : Numerical simulation results..

Figure 17 :

 17 Figure 17: Collision frequency with respect to the bed solid volume fraction (left), nondimensional collision frequency (right) f * coll = f coll

  .

	Liquid phase	ρ f	1400 kg/m 3
		µ f	3.8 × 10 -3 P a.s
	Fluidization velocity	U F	0.17/0.15/0.12/0.09/0.073 m/s
	Particles	ρ p	2230 kg/m 3
		d p	6 × 10 -3 m
		V t	0.24 m/s
		Re t	530
		St t	5.3
	Fluidization law

Table 1 :

 1 Phase properties and fluidization parameters

	3. Numerical model
	Details of numerical approach and validation test cases are given in Vincent

Acknowledgements

This work was granted access to the HPC resources of CINES under the allocation x20142b6115 made by GENCI (Grand Equipement National de Calcul Intensif) and of CALMIP under the allocation p0633 (Calcul en Midi-Pyréenés).

This work has been supported by the research federation FERMaT (FR CNRS 3089).

Appendix A. Definition of the statistical operators.

The average operators of a quantity φ are defined in a Eulerian or Lagrangian way, whether φ describes the behavior of the continuous or discrete phase. The different average operators are explicitly defined in this section.

-Lagrangian averages-

• Arithmetic average on the ensemble of particles

where φ n is a variable associated to the n th particle. N p is the total number of particles, and ∆T is the total simulation time. The time step for statistical calculations of the Eulerian phase is 50 times the simulation time step.

• Average over cylindrical shells of height h b and thickness ∆r such that

where φ i is the value of the variable φ of the continuous phase defined on the i th cell.

• Average over a disck of diameter D and thicknes ∆z,

where x i is the position of the i th cell.

• Phase average over the Eulerian grid

where χ i is the solid volume fraction of the i th cell.

• Phase average on a plan P = {x, x • e y = 0}.

[φ] plan =


Jean-Franc ¸ois Patri
email: jean-francois.patri@gipsa-lab.fr

Pascal Perrier

Julien Diard

Bayesian modeling in speech motor control: a principled structure for the integration of various constraints

Keywords: Speech motor control, Bayesian modeling, Hypo/Hyperspeech, Dynamical constraints

Speaking involves sequences of linguistic units that can be produced under different sets of control strategies. For instance, a given phoneme can be achieved with different acoustic properties, and a sequence of phonemes can be performed at different speech rates and with different prosodies. How does the Central Nervous System select a specific control strategy among all the available ones? In a previously published article we proposed a Bayesian model that addressed this question with respect to the multiplicity of acoustic realizations of a sequence of phonemes. One of the strengths of Bayesian modeling is that it is well adapted to the combination of multiple constraints. In the present paper we illustrate this feature by defining an extension of our previous model that includes force constraints related to the level of effort for the production of phoneme sequences, as it could be the case in clear versus casual speech. The integration of this additional constraint is used to model the control of articulation clarity. Pertinence of the results is illustrated by controlling a biomechanical model of the vocal tract for speech production.

Introduction

A number of various experimental studies have suggested that speech is a skilled serial-order motor task that is adapted to achieve time series of goals within a timing that does not allow any on-line cortical processing of feedback signals. This is illustrated for the auditory feedback by the absence of impact on speech production of a feedback delayed by up to 75 ms [START_REF] Stuart | Effect of delayed auditory feedback on normal speakers at two speech rates[END_REF]. In addition, the speech motor system is highly redundant with many available degrees of freedom, which makes the inference of motor commands from physical signals an "ill-posed" inverse problem [START_REF] Atal | Inversion of articulatory-to-acoustic transformation in the vocal tract by a computer-sorting technique[END_REF]. These degrees of freedom are used in different ways to deal with sequence planning and variations in the conditions of articulation. For instance, clear anticipatory behaviors have been observed in articulatory and acoustic patterns associated with speech sequences [START_REF] Farnetani | Coarticulation and connected speech processes[END_REF].

To deal with this complexity, motor control models classically consider a feedforward control scheme for ongoing speech production under normal conditions. This is associated with a feedback controller enabling, with a certain delay, a correction of the motor commands in case of inadequate feedforward commands or in presence of external perturbations [START_REF] Wolpert | Internal models in the cerebellum[END_REF]. In this context, speech planning, which aims at solving the "ill-posed" inverse problem by finding the motor command patterns adapted to the production of a speech sequence, has been classically modeled within an optimal motor control framework. This approach has been shown to generate results in close agreement with experimental data, in particular in terms of adaptation to perturbations [START_REF] Guenther | A theoretical investigation of reference frames for the planning of speech movements[END_REF] or in terms of anticipatory behavior [START_REF] Guenther | Speech sound acquisition, coarticulation, and rate effects in a neural network model of speech production[END_REF][START_REF] Perrier | Speech planning for V1CV2 sequences: Influence of the planned sequence[END_REF].

However, criticisms of this approach include key issues in cognitive science, such as the nature of the neuro-physiological mechanisms likely to be associated with cost computation and cost minimization, or the inability to account for the wellknown token-to-token speech variability [START_REF] Perkell | Variability in production of the vowels /i/ and /a[END_REF]. This last drawback is inherent to the feedforward optimal control scheme, since this scheme basically cancels possible variations along the degrees of freedom directions, by specifying a unique optimal solution to the control problem.

In this context we have published an alternative approach by formulating feedforward optimal control in a Bayesian modeling framework [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF], which has been implemented for the control of a 2D biomechanical model of the tongue [START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF]. The methodology is based on Bayesian Programming [START_REF] Lebeltel | Bayesian robot programming[END_REF][START_REF] Bessière | Bayesian Programming[END_REF][START_REF] Bessière | Probabilistic Reasoning and Decision Making in Sensory-Motor Systems[END_REF], which proposes a structure for the construction of Bayesian models. Bayesian modeling, in a nutshell, solves inference problems by computing probability distributions for the solutions instead of proposing specific values; this allows to solve ill-posed problems while conserving token-to-token variability in a principled way [START_REF] Colas | Common Bayesian models for common cognitive issues[END_REF]. It also preserves the basic principles underlying the search for optimality without being explicitly driven by the minimization of a cost, relying instead on probabilistic evaluation of all possibilities.

The approach was illustrated by reformulating GEPPETO, an existing optimal control model for speech production planning developed in the lab [START_REF] Perrier | Modeling the production of VCV sequences via the inversion of a biomechanical model of the tongue[END_REF], into the Bayesian modeling framework. We have shown that these two models are nested, with optimal control as a special case of the Bayesian model: indeed, the Bayesian model is simplified into the optimal control model, when the inferred control commands are strictly limited to those with maximum posterior probability. Variability in the Bayesian model is formally generated by assuming that control is performed by sampling motor commands randomly according to the distribution that solves the inference problem.

An interesting aspect of the proposed formalism is its coherence for dealing with multiple constraints. This is illustrated in the present work by describing an extension of the Bayesian model that takes into account an important additional constraint associated with the trade-off between accuracy and effort at fast speaking rate, which has also been integrated in GEPPETO [START_REF] Winkler | A model of optimal speech production planning integrating dynamical constraints to achieve appropriate articulatory timing[END_REF]. While our recently published model specified the con-straints only in terms of auditory goals to be reached at the different phonemes of a speech sequence, the proposed extension integrates the specification of the level of effort (weak, medium, strong) involved in the production of the sequence. Depending on this effort level, different probability distributions are associated to the achievement of the task. The inferred command patterns will be assessed on the biomechanical model of the tongue with a special focus on how increasing the speaking rate affects the articulatory and acoustic accuracy in phoneme production.

The paper is divided into four sections. Section 2 summarizes the main hypotheses involved in the formulation of the Bayesian model, with an emphasis on the additional constraint that is to be included. From these ingredients the extension of the Bayesian model is introduced in Section 3. Section 4 presents the results of the inference and the assessment of the inferred command patterns with the biomechanical tongue model. Finally, the strengths of the proposed modeling framework are discussed with respect to its capacity to deal with different kinds of physical constraint applied to speech production.

Methodology

Basic control of the biomechanical model

This section summarizes all hypotheses considered in the model. The first half of these hypotheses are exactly the same as those used to define our recently published model [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF]. The second half concerns the introduction of the additional totallevel-of-force constraint.

As in GEPPETO, we are concerned with the selection of control variables for the production of sequences of phonemes. We only consider phonemes {/i/, /e/, /E/, /a/, /oe/, /O/, /k/} that do not require jaw or lip movements since the version of the biomechanical model that is used for simulations only includes an account of the tongue [START_REF] Payan | Synthesis of VV sequences with a 2D biomechanical tongue model controlled by the equilibrium point hypothesis[END_REF][START_REF] Perrier | Influences of tongue biomechanics on speech movements during the production of velar stop consonants: A modeling study[END_REF]. Control variables correspond to the six control parameters, called λ [START_REF] Feldman | Once more on the equilibrium-point hypothesis (λ model) for motor control[END_REF], used to pilot the biomechanical model. The first three formants of the acoustic signal are assumed to define the auditory space in which motor goals associated with phonemes are specified as convex target regions. Finally, the knowledge of the mapping between control variables and formant values is assumed to be stored in an internal model in the CNS [START_REF] Guenther | A theoretical investigation of reference frames for the planning of speech movements[END_REF]. In our study it is implemented by a radial basis function (RBF) network learned from more than 2.10 4 simulations covering the whole motor space [START_REF] Poggio | A theory of networks for approximation and learning[END_REF]. The internal model is considered to be "static" as it associates commands and outputs at targets. This first set of hypotheses was used to specify in our previous model the selection of control variables based on the auditory characterization of phonemes only. However, different configurations of λ control variables may result in the same tongue shape (and therefore the same acoustic signal), corresponding to an identical equilibrium configuration for different generated forces. The total level of force influences the capacity of the tongue to satisfy the speech requirements associated with increasing speaking rates. GEPPETO characterizes every articulatory configuration at the targets with a corresponding level of effort. These levels of effort are associated with the levels of total muscle force that are categorized in three levels, "Weak", "Medium" and "Strong" [START_REF] Winkler | A model of optimal speech production planning integrating dynamical constraints to achieve appropriate articulatory timing[END_REF]. Since muscle force capacity is highly muscle dependent [START_REF] Pruim | Forces acting on the mandible during bilateral static bite at different bite force levels[END_REF][START_REF] Brand | The sensitivity of muscle force predictions to changes in physiologic crosssectional area[END_REF] and since phonemes are differentiated by the patterns of recruited muscles [START_REF] Honda | Organization of tongue articulation for vowels[END_REF][START_REF] Baer | Electromyography of the tongue muscles during vowels in /OpVp/ environment[END_REF][START_REF] Buchaillard | A biomechanical model of cardinal vowel production: muscle activations and the impact of gravity on tongue positioning[END_REF][START_REF] Waltl | An EMG study of the German vowel system[END_REF], for a given total level of effort, the actual muscle force involved depends on the considered phoneme. Thus, effort is not simply linked with the absolute amplitude of muscle force, but with the amplitude of the force relatively to the maximal capacity of the involved muscles to produce force. This is a way to not penalize phonemes, such as /i/ or /u/, that requires the activation of intrinsically strong muscles, such as the Genioglossus or the Syloglossus, as compared to phonemes, such as /a/ that is associated with a weaker muscle, the Hyoglossus. Here again, the knowledge of the control-to-force mapping is assumed to be stored in the CNS in a "static" internal model, which results from a learning process that generalizes the relation between motor commands and generated forces from a number of examples. It is implemented by a second RBF network.

Bayesian model

Model definition

Variables Variables correspond to the formal representation of all the relevant quantities selected for the description of the system. In the Bayesian framework, these quantities correspond to probabilistic variables. These variables are extracted from the control scheme presented in Section 2:

Φ is a discrete variable representing phonemes. It corresponds to the set of phonemes specified in Section 2.

S is a continuous vector variable in the auditory space described as the 3-dimensional formant space: S = (F1, F2, F3).

M is a 6-dimensional continuous vector variable representing the motor commands controlling the tongue:

M = (λ1, • • • , λ6).
ν is a continuous scalar variable representing the amount of total muscle force generated.

N is a discrete variable representing the level of effort generated at the target articulation of a phoneme. Its elements are the 3 levels of effort considered in GEPPETO.

Decomposition Our aim is to completely specify the joint probability distribution P (M S Φ ν N ). To this aim, the joint probability distribution is first decomposed following the standard chain rule of probabilities. Then, assumptions related to dependencies linking different variables are exploited, in order to simplify the obtained decomposition. The resulting decomposition is:

P (M S Φ ν N ) = P (M ) P (S | M ) P (Φ | S) P (ν | M ) P (N | Φ ν). (1)
The assumptions leading to the terms in the first line of Eq. ( 1) are the same as for our previous model [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF]. Terms appearing in the second line of Eq. ( 1) correspond to the additional constraint included in the model. The simplifications leading to these terms are the assumptions that the total force level ν is uniquely specified by the motor variable M , and that the level of effort N is completely characterized by the knowledge of Φ and ν. It is worth noting that the global structure of the obtained decomposition results in the combination of the structure involved in our previous model with the structure relating muscle forces and related level of effort. Figure 1 illustrates the incorporation of these two structures into the current model. Parametric forms Once defined the decomposition for the joint distribution P (M S Φ ν N ) (Eq. ( 1)) has been defined, it is necessary to specify the form taken by each of its factors. P (M ), the probability of occurrence of a motor pattern M , is specified by a uniform distribution, since, in the absence of further knowledge, there is no reason to assume that some motor control patterns should be favored as compared to others.

P (S | M ) and P (ν | M ) represent probabilistic accounts of the knowledge of the motor-to-auditory and motor-to-force mappings respectively (see Section 2). Since these mappings uniquely specify the auditory and force outputs corresponding to motor variable M , they are described by Dirac delta functions located on the outputs of the corresponding RBF networks. These outputs are denoted by s * (M ) and ν * (M ) respectively.

P (Φ | S) corresponds to the probability of associating phoneme Φ with the auditory input S. It is specified in the same way as in our previous model [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF], by defining a sub-model that implements the inference of this categorization based on Gaussian descriptions of the P (S | Φ) distributions, describing the dispersion of auditory variable S for each phoneme.

P (N | ν Φ) characterizes the link between the level of effort N and the total level of force ν for a given phoneme Φ. It can be seen, for each phoneme, as the categorization of the generated forces ν into one of the three level of efforts N . As for P (Φ | S), this term is specified by a sub-model that performs this categorization based on the probability distributions P (ν | N Φ) of forces ν generated with respect to each level of effort N and for each phoneme Φ. P (ν | N Φ) are assumed to be Gaussian, with location parameters (means) specified by the effort constraint implemented in GEPPETO, and variances scaled by parameter κν , which controls the strength of the constraint.

Inference of motor commands M producing phoneme Φ with a desired level of effort N

Having specified the joint distribution P (M S Φ ν N ), we now formulate the question that is to be solved by the model. As the problem is to infer the control commands M producing a desired phoneme Φ under a desired level of effort N , the approach consists in using the model to infer the probability distribution over M , conditioned on the specified values of Φ and N . The corresponding distribution, P (M | Φ N ), is obtained by standard Bayesian inference and is given by:

P (M | Φ N ) = S,ν P (M S Φ ν N ) M,S,ν P (M S Φ ν N ) . ( 2
)
Since in Eq. ( 2) the denominator is constant for fixed N and Φ, we focus on its numerator and make the denominator implicit by using a proportionality symbol "∝". Using the decomposition given by Eq. ( 1) and including the constant P (M ) term in the proportionality symbol, we obtain:

P (M | Φ N ) ∝ S,ν P (S | M ) P (Φ | S) P (ν | M ) P (N | Φ ν) ∝ P (Φ | s * (M )) P (N | Φ ν * (M )), (3)
where the summation over S and ν is reduced to the terms ν = ν * (M ) and S = s * (M ) for which P (ν | M ) and 3) completely defines the distribution from which motor variables M can be selected in order to produce a specific phoneme with a desired level of effort.

Inference of motor commands producing a sequence of phonemes with a desired level of effort

Until now we have only considered the selection of motor commands for the production of isolated phonemes. The interest of taking into account the level of effort corresponding to the planned motor variable is made clear in the generation of sequences of phonemes. We previously proposed a Bayesian model planning motor commands for the production of sequences of phonemes, which only considered the auditory constraint [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF]. Based on the above methodology for the planning of a single phoneme, it is straightforward to include the additional effort constraint into the sequence planning model. Figure 4 represents the diagram corresponding to this extended sequence-planning model from which we can infer motor variables for the production of sequences of phonemes with different levels of effort. Variables are duplicated and indexed by their order in the intended sequence. The additional variable Cm implements a motor constraint that imposes the proximity of planned motor commands resulting in anticipatory coarticulation effects [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF].

Results

Inferred control patterns

The decision policy underlying the selection of a set of control variables for the production of phoneme Φ with effort level N consists in random sampling based on P (M | Φ N ). This sampling is approximated by the Metropolis-Hasting algorithm that performs a Markov Chain Monte Carlo (MCMC) random walk that converges to the desired distribution. Simulations were performed with Matlab's "mhsample" function with 20 chains of 2.10 4 samples each and a burning period of 10 3 samples. Tests revealed that further increasing the number of samples had no influence on the global shape of the obtained distributions.

In order to validate the performance of the model it is necessary to evaluate whether control samples M obtained from this sampling policy effectively result in intended auditory and force values that satisfy the constraints. Figure 2 gives an example of the results. It represents the histograms of the first three intended formant values corresponding to the samples M obtained for phoneme /i/ from P (M | Φ N ) at the three effort levels. These intended formant values where computed from the RBF network used for the internal models described in Section 2. It can be seen that the formants correctly distribute inside the target regions for all the categories of effort. This remains true for all other phonemes. Similarly, Figure 3 represents the histograms of the total intended muscle force obtained from the same set of samples M used in Figure 2. Again, these levels of intended total force were obtained from the RBF network. It can be seen that the forces correctly distribute according to the corresponding regions that define each level of effort. In summary, the Bayesian model correctly infers control samples M that jointly satisfy both the auditory and the effort constraints characterizing the speech planning task.

Sequence planning at various levels of effort

We consider for illustration the sequence /aiE/ planned with the N="Weak" and N="Strong" levels of effort. The set of control variables having the highest inferred probabilities are selected and the resulting tongue trajectories are generated by the biomechanical model presented in Section 2. Two speaking rates are implemented for each set of control variables by specifying a slow and fast transition rate between the motor commands of the first and second phoneme in the sequence, together with a long and a short duration of the second phoneme. Results are illustrated as formant trajectories in auditory space shown in Figure5. It can be seen that for a fast speaking rate, motor commands planned with the N="Weak" level of effort result in a formant trajectory that misses the auditory target region for the middle phoneme /i/ (Top left image). This is not the case for control variables planned with the N="Strong" levels of effort (Top right image). In a slow speaking rate however, these same two sets of control variables result in formant trajectories that both reach auditory target regions (Bottom images).

Conclusion

We have presented an extension of a recently published Bayesian model of speech motor control that plans control variables for the production of phoneme sequences [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF]. This extension shows how an additional constraint, related to the level of total muscle force generated in the production of phonemes, can be included in a modular manner into the model. We have shown that the model correctly infers motor commands that jointly satisfy both constraints characterizing the speech task: intended auditory outputs lie inside the desired phonemic regions and total level of muscle force satisfies the desired effort interval. The model was constructed as a reformulation of GEPPETO, an existing optimal control model that rests on Figure 5: Formant trajectories corresponding to the motor commands planned for the sequence /aiE/ with the N="Weak" (Left) and N="Strong" (Right) effort levels, produced with fast (Top) and slow (Bottom) speech rates. For a fast speaking rate, the trajectory planned with the N="Weak" effort level misses the auditory target for the middle phoneme /i/ (Top left image), but not for the N="Strong" effort level. Both commands produce trajectories that reach the auditory targets for slow speech rate (Bottom images). Red crosses indicate the intended auditory output of the motor commands.

the minimization of a cost function. We have shown in a previous paper [START_REF] Patri | Optimal speech motor control and token-to-token variability: a Bayesian modeling approach[END_REF] that the Bayesian model includes the optimal control approach of GEPPETO as the particular case of finding motor commands with the maximum posterior probability.

In this sense the Bayesian framework offers a more general approach that has valuable properties. For instance, the Bayesian approach results in a distribution of solutions while the optimal control approach gives a unique point-like value. These distributions allow us to better explore the space of solutions and study the possibility that patterns of variability may be attributed to the distribution of possible equivalent motor solutions.

We see that the proposed Bayesian Programming approach offers an interesting framework for the integration of multiple constraints. The flexibility and coherence of the Bayesian methodology allows us to further include additional constraints considered in speech motor control. In this context, current work is focused in the combination of constraints related to proprioceptive characterization of phonemes along with the present auditory representation.

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Communitys Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 339152, "Speech Unit(e)s", PI: Jean-Luc Schwartz).

Figure 1 :

 1 Figure 1: Incorporation of the additional effort-constraint into our previous Bayesian model [9].

Figure 2 :

 2 Figure 2: Diagram of the sequence-planning model including the effort constraint.

Figure 3 :

 3 Figure 3: Histograms of the total level of force corresponding to samples M obtained by the the probability distribution P (M |ΦN ) for phoneme /i/ and effort levels N="Weak" (Left), N="Medium" (Middle) and N="Strong" (Right). The vertical dotted lines indicates the borders of the regions characterizing the effort constraint of GEPPETO.

Figure 4 :

 4 Figure 4: Histograms of auditory consequences of samples M sampled from the probability distribution P (M | Φ N ) for phoneme /i/ and effort levels N="Weak" (Bottom), N="Medium" (Middle) and N="Strong" (Top). The vertical dotted lines indicate borders of the auditory regions characterizing the auditory constraint for this phoneme.