François Kneib
email: francois.kneib@irstea.fr

Thierry Faug
email: thierry.faug@irstea.fr

Gilles Nicolet
email: gilles.nicolet@irstea.fr

Nicolas Eckert
email: nicolas.eckert@irstea.fr

Mohamed Naaim
email: mohamed.naaim@irstea.fr

Frédéric Dufour
email: frederic.dufour@3sr-grenoble.fr

Force fluctuations on a wall in interaction with a granular lid-driven cavity flow

The force fluctuations experienced by a boundary wall subjected to a lid-driven cavity flow are investigated by means of numerical simulations based on the discrete element method. The timeaveraged dynamics inside the cavity volume and the resulting steady force on the wall are governed by the boundary macroscopic inertial number, the latter being derived from the shearing velocity and the confinement pressure imposed at the top. The force fluctuations are quantified through measuring both the autocorrelation of force time-series and the distributions of grain-wall forces, at distinct spatial scales from particle-scale to wall-scale. A key result is that the grain-wall force distributions are entirely driven by the boundary macroscopic inertial number, whatever the spatial scale considered. In particular, when the wall-scale is considered, the distributions are found to evolve from nearly exponential to nearly Gaussian distributions by decreasing the macroscopic inertial number. The transition from quasistatic to dense inertial flow is well identified through remarkable changes in the shapes of the distributions of grain-wall forces, accompagnied by a loss of system memory in terms of the mesoscale force transmitted toward the wall.

I. INTRODUCTION

Granular materials are ubiquitous in nature and daily life situations. Depending on the stress level applied, dense packings of grains can exhibit solidlike or fluidlike behaviour. The rheology of dense granular flows is a key question among others in granular physics that has attracted increasing attention in the last twenty years [1][2][3]. The still unsolved questions, in particular related to the transition between quasistatic and dense inertial regimes, are relevant to a number of physical problems. Two examples are the modeling of the processes at stake when a full-scale granular flow-such as an avalanche or a landslide-initiates in the release area or ceases in the runout zone on the one hand, and the optimal design of silos to guarantee an efficient transport of particles in food processing or mining industry on the other hand.

Though outstanding progress was recently made regarding the rheology of dense granular flows, most of the existing models were developed to predict the average flow with only limited attention paid to the fluctuating part (fluctuating trajectories of grains, velocity fluctuations, stress fluctuations) of the flow [3]. In a similar manner, many studies about the force experienced by objects immersed in dense granular flows [4] focused on the average force signal but only a few of them (see for instance [5]) addressed the problem of the fluctuating part of the force signal.

There exists a large body of statistical studies which are focused on spatial force variability, more particularly on the analysis of the distributions of forces, in static (unsheared) or slowly sheared granular media, thanks to laboratory tests [6][7][8][9][10][11], numerical simulations [12][13][14] and theory [15][16][17][18][19][20][21]. In contrast, there are only few studies that have tackled the problem of temporal and/or spatiotemporal fluctuations (the distinction between both being not so clear in a number of studies) in slow to fast flows of granular materials [5,14,[22][23][24][25][26][27]. A brief overview of those studies is given in appendix A.

The present paper proposes to investigate in detail the fluctuating part of the force experienced by a boundary wall subjected to a granular flow. Our study is based on an original system, namely the granular lid-driven cavity, which is simulated by the discrete element method (DEM). This follows a preliminary analysis of the timeaveraged dynamics of this granular cavity system that was presented in [28]. The force fluctuations experienced by the boundary wall are analysed with the help of a systematic characterization of the probability distributions of grain-wall forces, under different confinement pressures and shear velocities at the top of the cavity. Moreover, various spatial scales are analysed including particle microscale, wall macroscale, and mesoscale (intermediate between the two previous scales). The original system considered in the present paper allows to investigate, for the first time, how the distributions of grain-wall contact forces evolve over a wide range of both shearing velocity and confinement pressure imparted to a complicated dense granular-flow geometry, and considering different spatial scales from a single grain to the entire boundary wall.

The present paper is organized as follows. Following our initial numerical study on the time-averaged dynamics of the granular lid-driven cavity system [28], Sec. II gives a summary of the macroscopic boundary conditions and the microscopic parameters which we are using for this new numerical study. Section III recalls briefly the main results regarding the time-averaged dynamics and enriches them with analysis of the horizontal position of the vortex formed within the cavity, as well as of the local rheology in the granular bulk. Section IV tackles the temporal fluctuations of force by analysing in detail the autocorrelation of the force time-series on both the entire height of the wall and portions of the wall. Section V is devoted to the analysis of the grain-wall force distributions at micro-scale (particle scale), macro-scale (entire height of the wall) and meso-scale (portion of wall). Finally, the paper ends by discussing the main results and some potential implications for basic aspects of granular physics.

II. THE GRANULAR LID-DRIVEN CAVITY SIMULATED BY DEM

A planar assembly of spheres is trapped into a box made of four walls, as sketched in Fig. 1. The bottom horizontal wall is fixed and rough, while both lateral walls are fixed, smooth and spaced by a distance L = 5H. The upper horizontal wall, namely the lid, is rough and has an infinite length in order to allow for a uniform horizontal displacement. This choice was initially motivated by practical applications in geophysics which involve (i) a rough bottom, (ii) moving grains that entrain static grains (rough top wall), and (iii) obstacles with a smooth surface (lateral smooth walls). The top wall produces a constant shearing to the sample at velocity U through the x-axis direction, and remains horizontal as it is not allowed to rotate. A constant vertical confinement force F c is added continuously to the lid whose vertical position is computed through the DEM algorithm. Although the grains reaction to this force is not homogeneous along the lid-but rather exhibits an exponential shape (see figure 4 of our preliminary study on this system [28]), an equivalent macroscopic confinement pressure P can be defined from the system dimensions: P = F c /Ld, where d is the particle diameter. As the grains are forced to stay inside the cavity and the force profile is heterogeneous, a perpetual circulation takes place within the whole volume of the cavity, as further discussed in Sec. III of the present paper.

The following dimensionless macroscopic parameters are used to run the numerical simulations over a broad range of confinement pressure P and shear velocity U , namely N P and N U :

N P = P Φρ p gH , (1)
N U = U √ gd . (2)
The parameter N P is the ratio of the pressure P to the typical hydrostatic pressure associated with height H of a system under gravity. The gravity acceleration g = 9.81 m s -2 is used for convenience in order to facilitate the parallel with any potential real laboratory tests in the future, though the numerical system considered here is gravity-free. A constant macroscopic volume fraction Φ = 0.6 is considered here, which corresponds roughly to the random close packing of a three dimensional assembly of spheres of width d. The particle density was taken equal to ρ p = 2500kg m -3 . Our numerical simulations use spherical particles whose centers are forced to stay on a planar surface (plane (x, y) in Fig. 1). Any volume fraction measured in two dimensions (2D) is systematically transformed into an equivalent volume in three dimensions (3D), assuming Φ 3D = 2 3 Φ 2D if we compare a sphere of diameter d included in a cube of identical size d to a disk of diameter d included in a square of size d. The parameter N U is the ratio of the shear velocity U at the top to the typical velocity associated with one particle of diameter d under gravity. Those choices are also made for convenience to facilitate the parallel with any potential real laboratory tests in the future.

The cavity is initially filled of grains under gravity deposition in order to produce a dense granular packing of height H. Once the system reaches the static equilibrium, the gravity is set to zero and the shearing starts for a period of 10 seconds. As it will be defined in Sec. IV (see Fig. 7), the data recording starts one second after the shearing has started. One second typically corresponds to the maximal time needed for the system to reach a permanent regime in terms of force measured on the sidewall (see also more detail in [28]). In most of the simulations presented in the present study, H/d was taken equal to 30 with H = 0.3m and d = 1mm. Some results from numerical simulations with H = 0.3m but d = 0.3mm will be discussed in conclusion.

Following the analysis proposed by [1] at the local grain scale, if we consider the typical time t P = d ρ/P associated with the top confinement pressure P and the typical time t U = H/U equal to the inverse of the macroscopic shear rate, one can define the macroscopic inertial number I M :

I M = t P t U = d ρ/P H/U , (3)
where ρ = Φρ p holds for the density of the granular sample within the cavity.

In the present study, N P varied from 0.01 to 100 and N U from 1 to 20, which allows to investigate a wide range of granular flow regimes defined by I M ranging typically from 5 10 -4 to 0.6. A very slight variation of Φ * -the volume fraction actually measured in our simulationswith the boundary conditions in terms of U and P was observed. This point will be discussed in Sec. VI. In addition to I M , another dimensionless parameter that controls the system is the length L of the cavity relative to its height H. It has been shown in our preliminary study of this granular lid-driven cavity system [28] that the constraint L/H ≥ 5 should be respected to prevent an effect of the cavity length on the steady force experienced by the wall facing the shear displacement (namely the right sidewall shown in Fig. 1).

All the results shown in the present paper were obtained for L/H = 5, but their sensitivity to L/H was investigated by performing simulations with L/H = 7 and L/H = 10. This sensitivity analysis revealed that-as observed for the mean force, the force fluctuations remainded not influenced by L/H as soon as L/H ≥ 5.

The contact laws for normal (viscoelastic) and tangential (elastic with a Coulomb threshold) forces between particles and the corresponding micromechanical parameters (normal and tangential stiffnesses, damping coefficient accounting for inelasticity and related to restitution coefficient, interparticle friction), used in our DEM numerical simulations, are described in detail in appendix B (see also [28]).

III. TIME-AVERAGED DYNAMICS

The overarching goal of the present paper is to study in detail the force fluctuations experienced by the boundary wall that faces the direction of the shearing velocity at the top. In our previous study of this system [28], we reported an analysis of the results regarding the timeaveraged dynamics. Before going into the details of force fluctuations, this section proposes to recall some important results reported in [28] and to further extend some of the results concerning the time-averaged dynamics.

A. Steady vortex and mean force

The macroscopic inertial number I M , defined from the boundary condition [see Eq.( 3)], is the relevant parameter to quantify the flow inertial state and to predict the mean force experienced on the sidewall scaled by the force imposed at the top, namely F /(P Ld) (see figure 8 in [28]). More specifically, our previous study showed that the mean force on the right sidewall is entirely controlled by the changes in granular flow regimes, the latter regimes being governed by the macroscopic inertial number. In the quasistatic regime at low I M , F /(P Ld) is constant. In the dense inertial regime at intermediate I M , F /(P Ld) is a nearly-linear increasing function of I M . In the rapid regime at higher I M values, F /(P Ld) starts saturating. The measurement of time-averaged local granular flow velocity-fields revealed the formation of one single vortex occupying the whole volume of the cavity. This cavity-scale vortex is caused by the fact this system forces the grains to move within the cavity volume, without any possibility to escape it. Figure 2 shows examples of the streamlines inside the cavity for low (I M = 1.2 10 -3 ), intermediate (I M = 6 10 -2 ) and high (I M = 0.3) values of I M . The x-position of the vortex centre along the cavity was extracted from the streamlines (as shown by the vertical dotted lines in Fig. 2), thus allowing to deduce the horizontal distance between the vortex center and the right sidewall. The Fig. 3(a) shows scaled by the cavity length L as a function of I M . In the quasistatic regime (I M 10 -2 ), is relatively constant at an approximative value of L/2, meaning that the vortex is centered. Beyond I M 10 -2 , the symmetry of the velocity pattern is lost, as depicted by /L which starts decreasing sharply.

F/PLd I M <10 -2 I M >10 -2 c 1 -c 2 L , R 2 =0.96
FIG. 4. F /(P Ld) versus /L. The grey dashed line is a fit of Eq.( 5) for IM > 10 -2 (R 2 = 0.96 is the regression coefficient). This graph reveals that the mean force on the wall is entirely controlled by the distance of the vortex from the right sidewall.

The vortex centre is more and more shifted toward the right sidewall when I M is increased, thus producing an important horizontal asymmetry at the highest I M . Figure 3(b) depicts how the mean force on the wall F relative to the typical pressure force P Ld imposed at the top evolves with the macroscopic inertial number (see [28] for more details). The variation of both F /(P Ld) and /L with I M reveals the same transition from quasistatic to dense inertial regime. While comparing the two curves, a significant coupling can be detected, which becomes obvious when plotting F /(P Ld) as a function of /L in Fig. 4. F /(P Ld) and /L are linearly linked, meaning that the mean force on the wall scaled by P Ld may be deduced from the vortex position, and vice-versa. It is worth noting that the scattering at low F /(P Ld) in Fig. 4 is due to the very low velocity of the flow in this region, making difficult the identification of the vortex centre position.

In our previous study, we proposed an empirical scaling for F as a function of I M in the following form [28]:

F = P Ld   r 1 + (r 2 -r 1 ) 1 1 + I 0 M I M   , (4)
where r 1 , r 2 and I 0 M were constant fitting parameters. For the data shown in Fig. 3, a good fit is obtained with r 1 = 0.53, r 2 = 0.86 and I 0 M = 0.02. The interparticle friction does affect the values of r 1 , r 2 and I 0 , as it will be discussed in Sec. VI.

Considering that the mean force on the wall is controlled by the distance of the vortex center from the wall ( being a function of I M as displayed in Fig. 3(a)), we can specify another scaling that is given by the following relation:

F = P Ld c 1 -c 2 (I M ) L , (5)
where c 1 and c 2 are two constants that may depend on the micromechanical parameters of the grains. In this specific case c 1 = 0.88 and c 2 = 0.5, as shown by the dashed line in Fig. 4.

The above analysis showing the link between the vortex center position and the mean force on the right sidewall (that faces the shearing direction at the top) further extends the results reported in [28] concerning the mean dynamics of the granular lid-driven cavity. For given micromechanical grain properties, all macroscopic physical quantities measured in the cavity can be deduced from the macroscopic inertial number defined on the imposed shearing velocity and confinement pressure.

B. The local µ(I)-rheology

To further test the rheology, the local effective friction µ loc and the local inertial number I were measured within the whole volume of the cavity. Detailed steps for the calculation of those local quantities are provided in appendix C. Though the strain field is rather complicated inside the cavity (see the streamlines in Fig. 2), we observed that the strain and the stress tensors were generally well aligned within most of the cavity volume. Figure 5 shows an example for I M = 6.09 10 -2 . The patterns were similar at higher I M . At lower I M , the collinearity between stress and strain tensors was generally well verified, except close to the top and bottom boundary walls for the lowest I M . This point will need further investigation in a future work. Under the reasonable assumption of collinearity between stress and strain tensors, the local µ(I)-rheology, given by the relation [29], was analysed. It was found to be valid in most of the cavity volume whatever I M . The values found for the parameters of the µ loc (I)-friction law were µ 1 = 0.12, µ 2 = 0.40 and I 0 = 0.13 (see much more details in appendix C C). The difference between the µ loc actually measured and µ th was calculated. Figure 6 displays the results in terms of maps of (µ loc -µ th )/µ th , for two distinct yet very close values of I M ∼ 10 -2 and I M ∼ 6 10 -2 . The latter value corresponds to the first point shown in grey color in Fig. 3 and the former value holds for the last point in black color in Fig. 3.

µ th = µ 1 + (µ 2 -µ 1 )(1 + I/I 0 ) -1
At I M ∼ 10 -2 , a symmetrical bow-like pattern appears across the entire length of the cavity (see (a)). Similar symmetrical bow-like patterns (not shown here) were systematically obtained for I M below 10 -2 . Above the bow-like pattern and below it (apart from the two regions close to the two lateral walls), µ loc -µ th is nil, which means that the spatial region is well governed by the µ(I)-inertial rheology. For higher values of I M (see Fig. 6(b)), the symmetry and the bow-like shape were broken. In other words, the region inside the cavity where µ loc -µ th was nil could extend over the whole height of the cavity.

Note that a great number of values of µ loc were found to be lower than µ 1 for I M below 10 -2 (see Fig. 19 in appendix). This corresponds to the region in the center of the cavity in Fig. 6(a), forming the bow-like pattern across the entire length of the cavity, where (µ loc -µ th )/µ th is negative. These negative values of µ loc -µ th suggest that the granular bulk in this region is able to flow below the friction threshold µ 1 . Such a situation was highlighted in a couple of granular systems and explained by non-local effects [30,31]. The conditions to have this bow-like pattern over the entire length of the cavity and its link to non-locality in dense granular flows (see [32] and references therein for an overview of non-local models for granular flows) will need further investigation in the future. However, it can be said that such a symmetrical bow-like appears when the boundary macroscopic inertial number is low, thus showing the limits of the local µ(I)-rheology in the quasi-static regime. In conclusion, the crossover from a bow-like pattern extending across the entire length of the cavity (at low I M ) to a smaller pattern in the confines of the boundary wall (at large I M ) appears to show proof of the quasistatic to dense inertial transition in the system studied here.

Moreover, one can note that triangular regions close to the lateral walls are well-identified in Fig. 6. Note that the collinearity between stress and strain tensors is well verified inside those zones, also identified in Fig. 5. In those regions, µ loc -µ th is positive, thus meaning that local friction is higher than the value predicted by the local µ(I)-rheology. This suggests that the µ(I)-rheology is not a sufficient tool to estimate the magnitude of the local force at the walls. However, we will show in the rest of the paper that the macroscopic inertial number is a good indicator for the transition from quasi-static through dense inertial to collisional regimes in terms of force fluctuations.

IV. FORCE DATA AUTOCORRELATIONS

A. Data recordings

While analysing temporal fluctuations, data acquisition frequency requires a particular attention: it has to be large enough to be able to capture short-life force events, and simultaneously small enough to avoid practical memory allocation issues while recording trends over large durations. The DEM algorithm is not able to model oscillating phenomena that occur at frequencies higher than the typical natural contact frequency f c which is the inverse of the period of the damped oscillator corresponding to a contact between two grains. In the present statistical study the record frequency f rec was set to 10 kHz. Using the formulation described in [33] for the collision time as a function of the grain properties gives a f c ranging from 4.8 kHz to 320 kHz. The frequency corresponding to the DEM timestep was systematically set to 10f c , avoiding any instability issues. The overlap between f rec and the lowest values of f c may cause the recording of unintended effects caused by particle oscillations. It has been verified that neither the autocorrelations nor the distributions of the force time-series presented in the following statistical study are sensitive to f rec by comparing to simulations with a doubled (20 kHz) and a halved (5 kHz) record frequency.

B. Force time-series autocorrelation

The force time-series on the right sidewall of the cavity show high temporal fluctuations, as displayed in Fig. 7. This section aims at studying these fluctuations by analysing the autocorrelation of force time-series. The first information that will be extracted is the force periodicity. The second information arising from force autocorrelation is the critical autocorrelation time ∆t c , which can be interpreted as the time after which the system has forgotten a certain force value on the wall. In other words, it reflects the typical time during which a future state of the system keeps history of its past state regarding the force transmitted to the wall.

We denote by C(F (t), F (t+∆t)) the temporal autocorrelation function, where ∆t is the lag between two system states. The critical autocorrelation time ∆t c is defined with an arbitrary low autocorrelation threshold:

C (F (t), F (t + ∆t c )) = 0.15. ( 6
)
It is worthwhile to note that the above threshold value does not influence the conclusions presented here as they are qualitative, focused on main trends. Since the local autocorrelation function C(F (t), F (t + ∆t)) is likely to be a non monotonic function, it is expected to cross the threshold many times. In this case, because we focus on the initial autocorrelation decrease corresponding to memory loss and we want to avoid any effect of the signal noise, the lowest value of ∆t c is kept.

In the following, we present the autocorrelation of force time-series at macroscopic scale, i.e. on the entire height of the wall (Sec. IV B 1). Then the autocorrelation at mesoscopic scale is studied by splitting the wall height into ten portions of a few grains each (Sec. IV B 2). Because of the strong stress localisation taking place at the top right corner of the cavity (see more detail on this wedge effect in [28]), the shearing frequency at the roughness (particle) scale is transmitted to the right sidewall, as long as I M is not too high (typically smaller than 10 -2 ). Sinusoid-like oscillations are therefore observed, with a decreasing period while increasing N U . The inset of Fig. 8 shows a collapse of sinusoid-like oscillations when the autocorrelation function is plotted versus (U/d)∆t, the wavelength associated with these oscillations being d (the diameter of the grains composing the top wall roughness and shearing the sample).

For the highest I M , this effect is absent, which is in agreement with a more inertial regime for which shorttime collisions occur in addition to enduring frictional contacts able to transmit the shear force. The mean free path of grains increases, which makes impossible any continuous transmission of enduring contacts between grains. Again, we detect the transition from the quasistatic to the dense inertial granular regime occurring around I M = 10 -2 , which is a value similar to the one extracted from the time-averaged dynamics in Sec. III.

The general shapes of the autocorrelation functions shown in Fig. 8 resemble a lot the ones found by Geng and Behringer [5] in their experimental study of an intruder slowly dragged into a granular medium (see figure 8 in [5]). We generally observe that C(F (t), F (t + ∆t)) drops quickly (exponentially) to zero over a time scale ∆t c and then fluctuate around zero. These fluctuations are well explained in our system by the typical frequency U/d associated with the shearing velocity at the top and the grain size (see discussion above). We were not able to find a clear dependency of the critical time ∆t c on either the macroscopic inertial number I M or any input parameter such as the shear velocity of the lid, or the confinement pressure at the top. In the following, we focus on portions of the wall. C(F(t),F(t +∆t))

I M 0.0012|N U =1.0|N P =25.0 I M 0.0061|N U =1.0|N P =1.0 I M 0.0122|N U =10.0|N P =25.0 I M 0.0913|N U =15.0|N P =1.0 I M 0.3043|N U =10.0|N P =0.

Autocorrelation for portions of the wall

The force time-series on the right wall is the result of the cumulative contact forces applied over the entire height of the wall. Consequently, the strong spatial heterogeneity of the granular lid-driven cavity system may require a more localized analysis by taking into account the position at which individual forces are applied to the wall. We propose here to investigate the autocorrelation of force time-series on some portions of the wall. The entire wall height is split into ten slices i ∈ {1...10} of identical size. The index i = 1 refers to the bottom slice and i = 10 refers to the top slice. Each slice has its associated force time-series from which the autocorrelation C i (F (t), F (t + ∆t)) is then computed.

Figure 9 shows the autocorrelation of force time-series for each of the ten identical wall portions, for three values of the macroscopic inertial number. The signal periodicity already discussed above is still observed at the frequency U/d-associated with the roughness of the top wall moving at velocity U . This is particularly clear in the simulation at I M = 1.2 10 -2 for which each temporal autocorrelation curve shows this period, whatever the vertical slice i considered along the wall. This indicates that the periodic fluctuations at the top can be transmitted to the lowest position at the bottom of the cavity. This result is consistent with the fact that we could identify a small, yet non-zero, spatial autocorrelation at I M = 1.2 10 -2 (not shown here). At the lowest I M = 1.2 10 -3 , the periodicity caused by the top roughness is present at the highest portion of the wall (i = 10) but quickly fades while going deeper into the sample along the wall, and even disappears for the lowest position at the bottom of the cavity. In contrast to the intermediate I M , the periodic fluctuations from the top wall cannot be transmitted over the entire depth of The full circles show <∆tc > that is the mean over the ten ∆t i c values. The data were extracted from the autocorrelation functions plotted on Fig. 9 the wall. This observation is in accordance with the fact that at I M = 1.2 10 -3 we did not find any spatial autocorrelation (not shown here). There might exist some zones near the cavity top inside which the grains have enough time to rearrange locally, thus being able to relax the high stress caused by jamming close to the wedge. At higher I M , the periodic fluctuations associated with the roughness of the top wall are lost whatever the vertical position along the wall, thus further confirming the transition toward a much more inertial regime for which the increase of the mean free path between grains prevents the transmission of those fluctuations (see Sec. IV B 1).

Figure 10 shows the mean value of autocorrelation time, < ∆t c >, averaged over all portions of the wall. Though there is some data scattering regarding ∆t i c (see grey-colored crosses on Fig. 10), this plot demonstrates that there exists a I M below which <∆t c > is not zero and beyond which <∆t c > vanishes. In other words, it means that beyond a I M ( 10 -2 ) it is not possible to predict a future system state from the past one, thus indicating that the memory of the system is completely lost. This observation was still valid by increasing L/H from 5 to 10 (not shown here). The granular lid-driven cavity system suggests here a remarkable transition from the quasistatic to the dense inertial granular regime, which is characterized by a total loss of system memory (not intrinsic to the material) in terms of the force transmission from the top wall toward the right sidewall.

V. FORCE DISTRIBUTIONS

The distributions of force time-series on the sidewall give crucial information on fluctuations, as they quantify the probability of each force value the sidewall may expe-rience. In this section, the analysis of force distributions concerns three spatial scales: the micro-scale focused on individual grain-wall contact forces, the macro-scale focused on the total force on the wall, and the meso-scale focused on the force experienced by portions of wall. For each scale, the distribution response to the same wide range of I M as tackled in the previous sections is investigated. All distributions presented in this section were obtained from quantities that were scaled by their timeaveraged value.

A. Force distributions at micro-scale

The distributions of contact forces measured in our DEM simulations, for grains in contact with the right sidewall, characterize the typical forces repartition applied locally on the wall (grain scale). Recorded contact forces on the wall from all recording times are merged to form the sample, and then the probability distribution is computed. Let us note f an individual grain-wall contact force and f the overall mean grain-wall contact force of a simulation. In the following, f = f /f is the scaled grain-wall particle contact force.

Figure 11 shows examples of distributions of f obtained with three different values of I M but from several distinct dipoles (N U , N P ). The remarkable result is that the distributions collapse very well at each I M , revealing that the repartition of particle contact forces on the wall is entirely controlled by the macroscopic inertial number of the granular lid-driven cavity.

Figure 12 shows the distributions of f computed for several I M . We have tested a number of well-established PDFs to fit these curves, as well as theoretical or empirical PDFs discussed in the literature for granular contact forces, specifically the ones proposed by [6] [Eq. (9)] and [21] [Eq. ( 10)] (see appendix A). None of those PDFs was able to give conclusive results over the whole range of I M tested in our DEM simulations. The results were conclusive, over the whole range of I M , with a truncated log-normal distribution that reads as follows:

P( f ) = 1 S 1 f σ √ 2π( f + f 0 ) exp - [ln( f + f 0 ) -f µ ] 2 2 f σ 2 , (7)
where f > 0, f 0 > 0, and S is the normalization factor corresponding to the value at f 0 of the survival function of the (untruncated) log-normal PDF. The scale parameter f µ , the shape parameter f σ and the location parameter f 0 are monotonic functions of the macroscopic inertial number I M , as shown in Fig. 12(b).

In appendix D, we provide the data in a log-linear plot (Fig. 20(a)) in order to highlight that the results are reasonable overall but not perfect. In particular, the small forces are quite well reproduced whatever I M but the tails are not well captured at the highest values of I M . Density Apart from the fact that the best fits were obtained with the truncated log-normal distributions over the whole range of I M tested, it remains challenging to provide physically-based arguments to justify the use of such a PDF. A key result of the present study is that the parameters of Eq. ( 7) are found to be well-defined monotonic functions of I M (see Fig. 12(b)). This opens a path to predict the force distributions empirically, as long as the boundary macroscopic inertial number is known. A question then arises: does that key result too hold for the distributions of the total force on the entire height of the wall?

I M 0.0061|N U =1.0|N P =1.0 I M 0.0061|N U =10.0|N P =100.0 I M 0.0061|N U =5.0|N P =25.0 I M 0.0304|N U =5.0|N P =1.0 I M 0.0304|N U =10.0|N P =4.0 I M 0.0304|N U =1

B. Force distribution at macro-scale

Figure 13 shows the distribution of F , the total force time-series exerted on the sidewall scaled by its mean F , for seven values of I M ranging from 1.2 10 -3 to 0.6. At the lowest I M in the quasistatic regime, the distribution are nearly Gaussian whereas at the highest I M , the PDF has a nearly exponential shape. This result is consistent with the crossover from exponential to Gaussian PDFs generally observed at high confinement pressures in some past studies [34][35][36] (see also the brief review proposed in appendix A).

At intermediate I M in the dense inertial regime, the distribution has a more complex shape: at a first glance, it may appear as a distribution which would combine some properties of a Gaussian-like distribution and an exponential-like distribution. However, it was impossible to find a way to model the distributions of total force by an analytical form that combines both Gaussian and exponential distributions, over the entire range of I M tested here. The only PDF that could be fitted to the distributions of the total force on the wall is again the trun- cated log-normal distribution [see Eq.( 7)], as shown by the dashed lines drawn in Fig. 13.

In the appendix D, we provide the same data in loglinear plots (see Fig. 20(b)) to highlight the robustness of the fits obtained with the truncated log-normal function over a wide range of the boundary inertial macroscopic number. As for the distributions of f discussed in Sec. V A, it would be difficult to give a physical interpretation for the truncated log-normal PDF for the distributions of F . However, it is possible to analyse the evolution of the corresponding parameters with the inertial number. F µ , F σ and F 0 (the scale, shape and local parameters, respectively) are presented on Fig. 13(b), as a function of I M . These parameters follow monotonous paths with I M , allowing to predict the distribution of the force on the sidewall from the inertial number. This result was further confirmed (with the same set of simulations as in Fig. 11) because several distinct dipoles (N U , N P ) that give the same inertial number systematically led to the same distribution (curves not presented here). The Fig. 13(b) displays a jump in the parameters of the truncated log-normal distribution between I M = 10 -2 and I M = 6 10 -2 . Further looking at the curves Fig. 13(a), this jump occurs when the nearly Gaussian shape almost disappears due to the competition with the exponential shape at low forces, thus displaying an homogeneous repartition (kind of plateau) from F/F = 0 to F/F = 1. This significant change (for I M ∼ 10 -2 ) of the distribution of the total force on the wall, when the low forces and the mean forces have almost the same probability of occurrence, appears to be an additional marker of the quasistatic to dense inertial transition, concomitant with the total loss of system memory discussed in Sec. IV B 2.

C. Force distribution at meso-scale

Because of the heterogeneity of the granular sample induced by the cavity boundary conditions, it is interesting to investigate the force distributions at a mesoscopic scale-smaller than the cavity height and greater than the grain scale, as already done for the force autocorrelation (see Sec. IV B 2). We note F i the force timeseries on a portion of wall located at a given i, and F i its mean. The Fig. 14(a) displays the distributions of F i /F i obtained on ten wall portions of identical size, for I M 10 -2 . The distributions corresponding to the portions located in the center of the right sidewall collapse into one single curve, while the distributions measured close to the top and bottom right corners of the cavity differ from this master curve. This result reflects the spatial heterogeneity of the cavity over its height, in connection to three zones: a highly sheared zone of small height at the top, a small dead zone at the bottom and (in between) a flowing zone of much larger height in the center. An identical behaviour (curves not shown here) was observed for all values of I M . These three zones, extracted from the distributions of F i /F i , are fully consistent with the vertical profile of the time-averaged force measured on the wall, as depicted in figure 5 of our previous study [28]. The time-averaged force was rather homogeneous over the same (large) central zone of the wall. It is worth noting that the sizes of the bottom and top zones may be sensitive to the grain diameter relative to the wall height. Further simulations with different grain diameters would be necessary to study the potential influence on the size of the boundary (top and bottom) zones identified here.

In the following, the collapse of the F i /F i distributions, excluding the extreme (top and bottom) zones of the wall, is exploited to analyse how the distribution at meso-scale evolves with the macroscopic inertial number. This mesoscale corresponds to a size H/10 = 3d.

The Fig. 14(b) depicts the mean distribution over seven values of i ∈ [3; 9] (corresponding to the central zone of the cavity excluding the top and bottom highly inhomogeneous layers), for a wide range of I M . As observed for the force distributions at micro-and macro-scale presented previously (see Sec. V A and V B, respectively), the curves depict shapes that resemble a truncated lognormal PDF. While comparing the distributions over the wide range of I M , they appear however to be much closer with each other, meaning they exhibit a weaker dependency toward I M . It has been verified (curves not shown here) that the distributions at macro-scale for the entire wall did remain unchanged if the data near the bottom and top wedges was not considered. The weaker dependency of the distributions to I M at meso-scale is thus striking, as it cannot be explained by a wedge effect only (this effect would be detected on the macro-scale otherwise). As a consequence, it can be concluded that similar meso-scale distributions can result in radically different macro-scale distributions through the summation relation that links the forces on wall portions to the force on the entire height of the wall. The latter observation may put emphasis on the crucial role of spatial dependencies between meso-scale force signals, that are controlled by the granular flow regime (through the value of I M ).

In particular, the Gaussian-like shapes shown on Fig. 13 obtained for quasistatic regimes (I M 10 -3 ) typically reflect the summation of independent meso-scale force signals. This is consistent with the fact that no spatial correlation was found at the lowest I M (as detected from the temporal autocorrelation functions at different i shown on Fig. 9(a)). Conversely, the evolution toward an exponential-like distribution at intermediate inertial numbers (I M 10 -2 ) on Fig. 13 reveals a spatial dependency between the meso-scale force signals. This is con-sistent with the spatial correlation detected between the temporal autocorrelation functions at different i shown in Fig. 9(b).

At I M 10 , the distribution displays a small amount of low forces and a nearly Gaussian repartition centered at about 0.7F i . Then, with the increase of I M , all distributions quickly converge to an exponential-like decrease with a high number of low forces. Finally, the force distribution on wall portions become independent of I M as soon as I M 10 -2 and i ∈ [3; 9]. This result allows once again the identification of a clear signature of the transition from quasistatic to dense inertial granular flow regime for the cavity system studied here.

VI. DISCUSSION AND CONCLUSION A. Time-averaged dynamics

Though there are still unresolved questions associated with complicated behaviours close to the lateral walls and with the formation of the bow-like patterns (probably caused by non-local effects), the analysis of the time-averaged dynamics in section III highlighted the role played by the µ(I)-rheology (proposed by [1,29]) in the granular force transmission toward a boundary wall. To further confirm it, we have analysed the slight variation of the height H of the cavity system. Up to now, a constant volume fraction Φ = 0.6 close to random close packing in three dimensions has been considered. The cavity system defined in our study is by construction a volume-free system, meaning that the volume fraction inside the cavity may vary. Indeed, though the variation of volume fraction of the granular bulk was small, one could measure it. The Fig. 15(a) displays how the volume fraction Φ * , actually measured in our DEM simulations and transformed into an equivalent volume fraction in three dimensions (see Sec. II), evolved with the boundary macroscopic inertial number I M . We observe a plateau (constant Φ * ) at the lowest I M , followed by a slow decrease of Φ * with log(I M ). This slight dilatancy of the granular bulk while increasing the inertial number is a robust result for a number of volume-free granular systems governed by the µ(I)-rheology, such as plane shear flows [37], free-surface granular flows down inclines [2], annular shear cell flows [38], etc. Finding this dynamic dilatancy law in the liddriven cavity system studied here further confirms the role played by the µ(I)-rheology. However, it should be kept in mind that the microrheology difference observed in some spatial regions of the cavity shows that the local µ(I)-rheology is not sufficient in particular close to the bottom and top boundary walls where non collinearity between strain and stress tensors was detected.

Note that the results in Fig. 15 the decrease of the grain diameter d can be explained by the boundary effects on the local volume fraction at both lateral and smooth walls of the cavity. This geometrical effect tends to decrease the volume fraction in the vicinity of smooth walls as the grains are forced to line up. A decrease in grain diameter fades the effect of the latter phenomenon on the macroscopic volume fraction. The Fig. 15(b) shows how the scaling between F /(P Ld) and I M , was influenced by the interparticle friction µ. The values of r 1 , r 2 and I 0 M used for the fitting function given by Eq. ( 4) are slightly changed when µ is divided by nearly two (see detail in Fig. 15(b)). Finally, the particle diameter d has a slight influence on the relation in the quasistatic regime, namely when I M is below 10 -2 . The influence of the grain properties (grain diameter, interparticle friction) will merit more attention in the future.

B. Grain-wall force distributions

Howell et al. [23] found exponential distributions at the smallest volume fraction while the distributions exhibited a peak at larger volume fraction (see appendix A). Considering the one-to-one relation between I M and Φ * shown in Fig. 15, the evolution of the distributions with I M depicted in Figs. 12 and 13 are in qualitative agreement with the evolution of the distributions with the volume fraction reported by Howell et al. [23] in their 2D granular Couette experiments.

While varying I M and the spatial scale of interest, the probability distributions of grain-wall forces which we found in the present study (summarized in Fig. 16) cover the wide range of shapes discussed in the literature overview proposed in appendix A. In their twodimensional granular Couette experiments, Howell et al. [23] compared the distributions at the particle scale to the ones averaged over a great number of particles, and found that at large volume fraction the distributions became nearly Gaussian. The latter observation appears to be consistent with the evolution of the distributions for I M = 1.2 10 -3 (largest volume fraction in this numerical study) when the spatial scale is increased, as shown in Fig. 16.

A key result of the present study is that the value of the boundary macroscopic inertial number and, to a lesser extent, the spatial scale (macro, micro, or meso) considered, are key inputs that contribute to predetermine the transition observed in the grain-wall force distributions. Note that doubling L/H (not shown here) did not change the distributions. All measured distributions could be modeled by truncated log-normal PDFs whose parame-ters were fully controlled by I M . The fits were very good at the wall-scale (see Fig. 20(b)) and reasonable at the particle scale (Fig. 20(a)).

C. Quasistatic to dense inertial regime

The analysis of the granular lid-driven cavity problem showed that the macroscopic inertial number I M fully controlled all the time-averaged quantities (the mean scaled force on the wall, the position of the vortex, and the bulk volume fraction) and the force fluctuations (through the analysis of probability distributions). Moreover, a clear transition was observed at I M around 10 -2 in the time-averaged dynamics of the cavity system. This transition was identified by an increase of the scaled steady force, a displacement of the vortex position toward the wall, a loss of symmetry in the |µ loc -µ th | maps, and a macroscopic dilation of the granular bulk. The autocorrelation of force time-series at meso-scale allowed to reveal a total loss of system memory (not intrinsic to the material) beyond a value of I M ∼ 10 -2 . Finally, we identified drastic changes in the evolution of the shape of probability distributions beyond a I M once again close to 10 -2 . All these results are interpreted as a salient signature of the transition from quasistatic to dense inertial flow regime in granular media. However, predicting the force fluctuations on the wall over a wide range of I M remains a challenging question, as evidenced by the complicated patterns formed close to the boundary walls (Fig. 6).

ACKNOWLEDGMENTS

This work has been partially supported by the LabEx Tec21 (Investissements d'Avenir: grant agreement No. ANR11-LABX-0030).

Thierry Faug and Mohamed Naaim are grateful to the financial support by the People Programme (Marie Curie Actions) of the EU 7th FP under REA Grant Agreement No. 622899 (FP7-PEOPLE-2013-IOF, GRAINPACT). François Kneib and Thierry Faug would like to thank Pierre Jop for instructive discussions. The authors thank two anonymous referees for their insightful comments and suggestions on this paper.

Appendix A FORCE DISTRIBUTIONS IN GRANULAR MEDIA

The present appendix provides an overview of the literature about contact force distributions in granular media for static (unsheared) or slowly sheared systems, as well as for fast flow systems.

For static granular packings, the q-model was early proposed and successful to reproduce the inhomogeneous interparticle contact force distributions observed in the pioneering experimental and numerical studies [15,39].

In particular, the q-model is able to predict the remarkable exponential decay at large forces generally measured in static granular matter. The high probability (compared to a Gaussian distribution) of having forces much larger than the mean is generally associated with the ability of granular materials to develop long chainlike structures, the so-called granular force chains, that can support the large forces [13]. The probability density function (PDF) of interparticle contact forces f predicted by the q-model in static granular media under gravity takes the following form [15,39]:

P(f /f ) = k k (k -1)! f f k-1 exp -k f f , ( 8
)
where f is the mean and k is the number of downward neighbour particles considered (see much more details in [15,16,39]). Considering k = 2-3 generally gives good results for the forces larger than the mean [12,34]. However, thanks to technological progress in the field of force sensors, a number of experimental studies identified a range of complicated shapes of contact force distributions when approaching forces much below the mean (see [12] for a detailed summary of those studies, before 2000, about distributions at weak forces). In particular, the presence of a plateau, followed by a slight increase at the smallest forces, was identified under certain circumstances. In order to fit that more complicated shape of contact force distributions, the empirical following functional form was proposed in [6]:

P(f /f ) = a 1 -b exp - f 2 f 2 exp - βf f , (9)
where a, b and β correspond to the PDF parameters fitted on the experimental measurements made by [6]. A slight modification of Eq. ( 9) was proposed by [8] to interpret their experimental data on the effect of both packing order (disordered packings versus highly orderedcrystalline-configurations) and interparticle friction on P(f ). Finally, it is worthwhile to note that log-normal distributions were reported in some studies [40,41].

As discussed in a recent review in [21] and a number of references therein, the key features of the contact force distributions in static or very slowly sheared granular packings can be summarized as follow: i) the distribution functions fall-off exponentially at large forces, ii) a small peak, or plateau, is observed below the mean force, and iii) the vanishingly small forces remain highly probable. In his review about granular force transmission in static granular packings, Radjai [21] proposed an elegant model for contact force distributions and derived an analytical expression for the density function able to predict the three aforementioned features:

P(f ) = β 0 (1 + γ 0 ) γ 0 exp(β 0 f ) [1 + γ 0 exp(β 0 f )] 2 , ( 10
)
where β 0 (homogeneous to a force) and γ 0 (dimensionless coefficient) are the PDF parameters. Under the normalization f = 1, β 0 and γ 0 are linked by the relation

β 0 = (1 + γ 0 ) ln[(1 + γ 0 )/γ 0 )].
Though further studies are needed to relate γ 0 to actual physical properties of the grains, varying γ 0 allows to cover a wide range of distributions with or without the presence of a peak [21].

Among the great number of experimental, numerical and theoretical results concerning the shape of force distributions in static packings or slowly sheared granular media, a key result is that the distribution at small forces (plateau versus peak, maximum value, non-zero value at vanishingly small force, etc.) is found to be very sensitive to the granular sample preparation and shear history the system experiences [12,21]. While interpreting distributions of contact forces in static granular packings, Antony [12] concluded that it is required to pay a due attention to the shear strain level and any other quantities related to shear history, such as the volume fraction. A theoretical study on static packings proposed by [19] predicted a broadening of the distribution while increasing the shear stress level, moving from distributions with a peak below the mean value to exponential distributions. Twodimensional granular packings under anisotropic stresses were studied theoretically by [20], showing that an increase in the stress anisotropy produced a transition from distributions with a peak below the mean to exponential distributions.

A few studies observed the possibility of a crossover in the shape of P(f ) from exponential to Gaussian at large confining pressure [34][35][36]. This change in the shape of P(f ) was attributed to very large deformations of particles [8]. The effect of particle stiffness on distributions of contact forces in static packings was investigated by Erikson et al. [9]. They found that beyond a high threshold-around 30% of deformation, the distributions became sensitive to the particle stiffness. Decreasing the particle stiffness led to an increase in the peak and of the slope of the exponential decay at large forces, thus producing a distribution likely to evolve toward a Gaussian distribution. The increase of the peak by decreasing the particle stiffness was also derived from the theory in [19] (note that this theory predicts an increase of the non-zero value at f = 0). In the limit of rigid grains (as addressed in the present study), one can conclude that the distributions are very weakly-not to say not-affected by the particle stiffness.

The pioneering laboratory measurements of contact forces distributions in granular media (see for instance [6-8, 15, 36]) were confined to contact forces between a grain packing and a boundary wall. It is worthwhile to note that some studies paid attention to the cross-comparison between the grain-grain (inside the bulk) contact force distributions and the grain-wall contact force distributions. No significant difference was found for static packings [18,19,42], suggesting that the results for the distributions at the boundary walls may generally be extended to the grain-grain contact forces.

For slowly sheared granular packings, the effect of a slight variation of the volume fraction on the shape of distributions was evidenced by the two-dimensional experiments on a granular Couette geometry [23]. The authors observed exponential distributions at the smallest volume fraction (around 0.78 in their 2D system), while the distributions exhibited a peak at larger volume fraction (around 0.81). Interestingly, Howell et al. [23] provided a comparison between the distributions at the particle scale to the distributions averaged over a great number of particles (∼ 260 in their study). For the latter, they found that at large volume fraction the distributions became nearly Gaussian. The transition from exponential distributions to narrower distributions with the presence of a peak while increasing packing density was also reported in a theoretical study based on the analogy with supercooled liquids and foams [43].

A study on relatively fast granular flows (a silo discharge) identified the effect of the velocity (flow rate at the exit of the silo) on the shape of distributions [24]: the authors observed a broadening of the distribution at large velocities. A similar conclusion was drawn from the experiments on objects dragged into a static granular medium in an annular cell by Geng and Behringer [5], who observed broader distributions at large rotation speeds of the annular cell. The theoretical study of O'Hern et al. [43] based on the analogy with supercooled liquids and foams predicted the broadening of distributions with the increase of the shear stress.

Appendix B MICROMECHANICAL PARAMETERS FOR DEM

The grain-grain and the grain-wall interactions are handled by a viscoelastic contact law for the normal force and an elastic force with a Coulomb threshold for the shear force. The normal contact force F n and the shear contact force F s are expressed as:

     F n = max(0, k n δ n + c n δn )n, dF s = (k s δs dt)s, |F s | ≤ µ|F n |, ( 11
)
where n is the normal of the contact plane, s is the unity vector along the shear direction (n • s = 0), k n and k s are the normal and tangential contact stiffnesses, δ n is the normal penetration depth, δs is the tangential displacement increment, µ is the local friction coefficient, c n is a damping factor that accounts for inelasticity of particles and is related to the restitution coefficient of particles, and dt is the timestep.

For a contact between two identical spheres of diameter d, Young's modulus E and Poisson's ratio ν, the contact stiffnesses can be computed as:

k n = 1 2 Ed, k s = νk n . (12)
The contact law is governed by four physical parameters: E, ν, c n and µ. It has been discussed in a number of studies (see [28,44] and references therein) that E can be reduced to decrease the total time of calculation without changing the numerical results, as long as we stay in the limit of rigid grains. Such a limit of rigid grains is verified if, for instance, we take N 0 = E/(2P ) = 1.5 10 4 , where P holds for the macroscopic confinement pressure and N 0 is a dimensionless number. According to Eq.( 12), the grains stiffness was not kept constant but set with respect to the macroscopic pressure applied P and the above criterion. As a result, the values of k n typically ranged between 6.6 10 1 N.m -1 and 6.6 10 5 N.m -1 over the whole set of simulations. A sensitivity analysis to grain stiffness was performed and revealed that-as long as N 0 ≥ 1.5 10 4 , all the results presented in this paper do not depend on this parameter. The Poisson ratio ν was taken equal to 0.3. c n is set in the same manner as in [28], with a restitution coefficient e = 0.5. Most of the simulations were conducted with µ = 0.5 but we will shortly discuss some results with µ = 0.27 in the conclusion of the paper (see Sec. VI). In order to minimize crystallization on the one side (low polydispersity) and to prevent migration and segregation mechanisms on the other side (high polydispersity), the grain diameters are homogeneously chosen between d(1 -d disp ) and d(1 + d disp ) where d disp = 0.15.

Appendix C FROM DISCRETE TO CONTINUUM MODELING A Kernel smoothing

As the granular sample is made of discrete particles involving discontinuities, the construction of classical fields often used in continuous mechanics require a particular attention. In this paper, a spacial kernel smoothing is used, as described in Fig. 17.

A 2D set of points m j spaced by d constitutes a regular mesh over the whole cavity. The physical quantities are computed at each mesh point from the particles quantities, according to the particle-point distance. Thus, any (scalar or tensor) quantity Q i which is associated with each individual sphere S i can also be associated with all individual mesh points m j , forming a continuous field. This quantity at the point m j , noted Q j , is the result of a spacial kernel smoothing with a Gaussian kernel:

Q j = 1 i N (D ij ) i Q i N (D ij ), (13)
where D ij are the distances between the spheres centers and the mesh points, and N is the Gaussian function of mean 0 and standard deviation d/2. It is worth noting that the use of a Gaussian kernel is arbitrary and any classical kernel such as a simple rectangular function should not modify the results. The choice of the standard deviation follows the result presented in [45]: the coarse-graining width should be of the order of magnitude of the grains diameter to avoid any oscillation effect and at the same time to limit sidewall effects.

B Local stress and strain rate tensors

The per-particle stress tensors and strain rate tensors used in the current paper are computed thanks to the algorithms that are part of the YADE-DEM open-source code (see [START_REF] Šmilauer | Yade Documentation[END_REF] and [START_REF] Catalano | [END_REF]). The space inside the cavity is discretized with the help of a Voronoï tesselation algorithm, as shown by the dotted lines in Fig. 18. In this way, it is possible to compute quantities based on an equivalent continuum of the cavity space. Each sphere is then associated with a list of neighbours forming a bounding polyhedral contour C with their centroid positions, and a bounding volume V σ which does not overlap with the neighbours.

Making the assumption that the spheres are at static equilibrium, the following Love-Weber expression can be applied to each sphere in order to compute their associated local stress tensors based on the contact forces:

σ = 1 V σ cn f c ⊗ l c , (14)
where c n is the ensemble of contacts on the sphere, f c are the contact forces and l c are the vectors linking the sphere center to the contact points. The local strain rate tensor D is defined as D = ∆t , where is the strain tensor computed from two system states shifted by a short lag ∆t. For a particular sphere, the displacements of all neighbours during ∆t are computed. Then the average displacement gradient < ∇d r > is obtained from the integration of the displacement d r along C (see [START_REF] Catalano | [END_REF] for more details):

< ∇d r >= 1 V D C d x ⊗ nds, (15)
where V D is the volume associated with the contour C , and d r is linearly interpolated on the segment between two successive neighbours. Finally, the strain tensor is given by the symmetric part of < ∇d r >.

C Inertial number and effective friction coefficient

The stress tensors σ and strain rate tensors D computed according to Sec. C B enable the computation of the spatial fields of inertial numbers I and effective friction coefficients µ loc . The following formulations are the same as described in [48] and are applied at each field point m.

The norm ||A|| and the deviatoric component A of a matrix A are defined as:

||A|| = Tr(AA T ) 2 , (16)
A = A - Tr(A) 3 I 3 , (17)
where I 3 is the identity matrix of size 3. The inertial number I and the effective friction coefficient µ loc are:

I = d √ ρ p ||D || √ p , (18)
µ loc = Tr(σ D ) 3

1 p||D || . (19)
In the above relations, p denotes the pressure that is defined from the decomposition of the Cauchy stress into the isotropic pressure and the deviatoric stress: σ = pI 3 + σ . Figure 19 displays µ loc as a function of I for each point m that belongs to a selected region of the cavity not too close to the boundary walls (see its definition in the caption of Fig. 19) and for all simulations of this study-over the whole range of I M tested. Finally, the empirical following form proposed by [29] is used to fit all the points, as drawn by the solid line on Fig. 19:

µ th (I) = µ 1 + (µ 2 -µ 1 ) 1 1 + I 0 /I , (20)
where µ 1 = 0.12, µ 2 = 0.40 and I 0 = 0.13. Note that a great number of points falls below µ 1 at low I M , suggesting that non-locality [30] is present. This point will need further investigation in the future.

FIG. 1 .

 1 FIG. 1. Geometry and boundary conditions for the lid-driven cavity simulated by DEM. The granular sample is trapped between four walls, the top one (the lid) applies the shear displacement U and the confinement pressure P (see text for details).

FIG. 2 .

 2 FIG. 2. Time-averaged streamlines within the cavity for three different values of IM . (a): IM = 1.2 10 -3 , NU = 1, NP = 25 ; (b): IM = 6 10 -2 , NU = 10, NP = 1 ; (c): IM = 0.3, NU = 10, NP = 0.04.

FIG. 3 .

 3 FIG. 3. (a): Distance of the vortex center from the right sidewall, , relative to the cavity length L, as a function of the macroscopic inertial number IM . (b): Mean force F on the wall relative to the pressure force P Ld imposed at the top, as a function of IM .

FIG. 5 .

 5 FIG. 5. Spatial maps of the principal strain (a) and stress (b) orientations averaged over time within the cavity volume: example for IM = 6.09 10 -2 , NU = 10, NP = 1.

FIG. 6 .

 6 FIG. 6. Maps of (µ loc -µ th )/µ th . (a): IM = 1.22 10 -2 , NU = 10, NP = 25 ; (b): IM = 6.09 10 -2 , NU = 10, NP = 1.

5 FIG. 7 .

 57 FIG. 7. Example of time-series of total force on the wall, obtained for IM = 1.2 10 -2 , NU = 10, NP = 25. The instantaneous force F is highly fluctuating over time. The horizontal dashed line depicts the value of the time-averaged force, named F , reached in the steady-state regime (after a short transient of about 1s). Inset: same data, zoomed in a time window of 0.05 s.

Figure 8

 8 Figure 8 depicts the autocorrelation of the force signal on the entire height of the right sidewall, for five values of the macroscopic inertial number I M .Because of the strong stress localisation taking place at the top right corner of the cavity (see more detail on this wedge effect in[28]), the shearing frequency at the roughness (particle) scale is transmitted to the right sidewall, as long as I M is not too high (typically smaller than 10 -2 ). Sinusoid-like oscillations are therefore observed, with a decreasing period while increasing N U . The inset of Fig.8shows a collapse of sinusoid-like oscillations when the autocorrelation function is plotted versus (U/d)∆t, the wavelength associated with these oscillations being d (the diameter of the grains composing the top wall roughness and shearing the sample).For the highest I M , this effect is absent, which is in agreement with a more inertial regime for which shorttime collisions occur in addition to enduring frictional contacts able to transmit the shear force. The mean free path of grains increases, which makes impossible any continuous transmission of enduring contacts between grains. Again, we detect the transition from the quasistatic to the dense inertial granular regime occurring around I M = 10 -2 , which is a value similar to the one extracted from the time-averaged dynamics in Sec. III.The general shapes of the autocorrelation functions shown in Fig.8resemble a lot the ones found by Geng and Behringer[5] in their experimental study of an intruder slowly dragged into a granular medium (see figure8in[5]). We generally observe that C(F (t), F (t + ∆t)) drops quickly (exponentially) to zero over a time scale ∆t c and then fluctuate around zero. These fluctuations are well explained in our system by the typical frequency U/d associated with the shearing velocity at the top and the grain size (see discussion above). We were not able to find a clear dependency of the critical time ∆t c on either the macroscopic inertial number I M or any input parameter such as the shear velocity of the lid, or the confinement pressure at the top. In the following, we focus on portions of the wall.

FIG. 9 .

 9 FIG. 9. Evolution of C i (F (t), F (t + ∆t)) over time measured on 10 portions i of wall, for three values of IM . In each graph, the horizontal dashed line shows the threshold of 0.15 that was used to derive ∆t i c . (a): IM = 1.2 10 -3 , NU = 1, NP = 25 ; (b): IM = 6 10 -2 , NU = 10, NP = 1 ; (c): IM = 0.3, NU = 10, NP = 0.04.

FIG. 10 .

 10 FIG. 10. Critical time ∆t ic (cross symbols) beyond which the force time-series becomes weakly correlated (C < 0.15)-for ten portions i ∈ [1; 10] of the wall height, as a function of IM . The full circles show <∆tc > that is the mean over the ten ∆t i c values. The data were extracted from the autocorrelation functions plotted on Fig.9

FIG. 12 .

 12 FIG.12. (a): probability distributions of contact forces for grains in contact with the wall, for different values of IM ; the dashed lines show a fit by a truncated log-normal distribution (maximizing the likelihood with the Nelder-Mead optimization method). (b): variation of the parameters of the log-normal distributions as functions of IM (see Eq.(7)).

FIG. 13 .

 13 FIG. 13. (a): PDF of the total force on the right sidewall for eight values of IM ranging from 1.2 10 -3 to 0.9. A truncated log-normal distribution can be fitted on the data (dashed lines). (b): the truncated log-normal parameters variation as a function of IM .

  FIG. 14. (a): PDF of force time-series on ten wall portions i for IM = 1.2 10 -2 , NU = 10, NP = 25 (each portion timeseries is scaled by its corresponding time-averaged force). The blue, orange and red curves respectively correspond to i ∈ [1; 2], i ∈ [3; 9] and i = 10. (b): mean distribution over i for i ∈ [3; 9] (yellow curves of (a)) for seven values of IM .

5 :FIG. 15 . 2 3

 5152 FIG. 15. Volume fraction Φ * (a) and scaled mean force on the sidewall F /(P Ld) (b) actually measured in the numerical simulations, as a function of IM : d = 0.001 and µ = 0.5 (blue curve), d = 0.0003 and µ = 0.5 (green curve), d = 0.001 and µ = 0.27 (red curve). (a): we show the equivalent volume fraction in three dimensions (3D) deduced from the volume fraction in two dimensions extracted from DEM (2D): Φ3D = 2 3 Φ2D (see discussion in Sec. II). (b): the fits proposed in [28] are drawn to show how µ affects the fitting parameters r1, r2 and I 0 M .

  FIG.16. PDFs of grain-wall forces at three different spatial scales for two extreme values of IM .

FIG. 17 .

 17 FIG. 17. Illustration of the kernel smoothing used for regular mesh fields calculations.

FIG. 18 .

 18 FIG.18. Illustration of the tesselation for the per-particle volume calculation.

FIG. 19 .

 19 FIG.19. Test of the local µ(I)-rheology for all simulations of this study, including all IM . The points come from a region of the cavity selected to approach the conditions of a simple shear flow. This region corresponds to the upper half region along y that is centered in a band of length L/2 along x.

  FIG. 20. Log-linear plots of the force distributions under different values of IM for the grain-wall contact forces at particle scale (a) (see detail in caption of Fig. 12(a)) ; for the mean force on the wall (b)(see detail in caption of Fig. 13(a)).

  .0|N P =0.04 I M 0.3043|N U =20.0|N P =0.16 I M 0.3043|N U =5.0|N P =0.01 I M 0.3043|N U =10.0|N P =0.04 FIG. 11. Probability distributions of contact forces for grains in contact with the wall. For each IM the curves are collapsing whatever the dipole (U, P ), thus demonstrating that the boundary macroscopic inertial number fully controls the grain-wall contacts distributions.

Appendix D FORCE DISTRIBUTIONS: LOG-LINEAR PLOTS

In this appendix, we provide the log-linear plots of force distributions for the contact force at particle scale (Fig. 20(a)) and the mean force on the wall (Fig. 20(b)), which correspond to the data already shown in linearlinear plots of Fig. 12(a) and 13, respectively.


J.-C David

S Chabod

D Doré

D Ene

B Rapp

D Ridikas

N Thiollière

In-Target Yields for Radioactive Ion Beam (RIB) Production with EURISOL

EURISOL DS (European Isotope Separation On-Line -Design Study) project is the European common effort in planning a next generation RIB factory able to deliver secondary beams up to 10 13 pps at energies up to 150 MeV u -1 . The proposed schematic layout of the facility is based on four target stations, three direct targets of 100kW of beam power and one multi-MW (MMW) target two stages assembly. Being produced via spallation the RIBs produced in the direct targets are mainly proton rich. While in the multi-MW target high intensity RIBs of neutron rich isotopes are produced by fission in actinide targets placed in the fast neutron spectrum given by a liquid metal spallation source. The purpose of this paper is to summarize a work carried out within Task 11 "Beam Intensity Calculations": estimation of the in-target yield intensities produced in the various target configurations. Benchmark studies were performed initially in order to verify the accurate description of the spallation models used by the MCNPX2.5.0 code and to choose the best options to be used for the present work requirements. Numerous calculations using MCNPX2.5.0 combined with the evolution code CINDER'90 were carried out to assess the performance of the direct targets. The production rates in the case of the MMW-fission targets were obtained with a given and fixed geometry (optimized to reach 10 15 fission s -1 ). Only fissile material, moderator and reflector were free.

Introduction

The European Community has launched the design study for a next generation RIB facility able to increase, by a few orders of magnitude, the exotic beam intensity and availability in Europe. In EURISOL [1], four target stations are foreseen, three direct targets of 100 kW of beam power and one multi-MW target assembly, all driven by a high-power particle accelerator (∼1 GeV). In this MMW station, high-intensity RIBs of neutron-rich isotopes will be obtained by inducing fission in several actinide targets surrounding a liquid metal spallation neutron source, whereas high energy spallation reactions will take place in the direct targets leading to proton-rich nuclei. A possible layout is given in FIG. 1.

Within this project a task were dedicated to the estimates of ion beam intensities. This paper reports the in-target yields obtained with several possible targets. It contains three main parts: benchmarking, direct target and fission target. The first part aims at defining which spallation models should be used and how well suited they are for EURISOL. In the second one numerous calculations were performed due to many possible configurations for the direct targets (material, dimensions, beam energy, etc.). The third one is devoted to the fission (or MMW) target. In this case the complicated geometry is fixed, but fissile material remains a free parameter, as moderator and reflector.

FIG.1. EURISOL facility layout (when EURISOL-DS project started in 2005)

Once obtained these yields should be used with extraction efficiencies to provide the realistic beam intensities.

Benchmarking

Although some valuable experience is already available from presently operating ISOL facilities world wide, the design of a new generation RIB factory requires specific and validated modelling tools. Then, the first step was to benchmark the model(s) possibly used on residue and neutron production with experimental data.

In MCNPX2.5.0 [START_REF]MCNPX User's Manual, Version 2.5.0[END_REF] we can use 10 spallation models (or to be more precise, model combinations, namely Intra-Nuclear Cascade model ⊕ Deexcitation model). We have begun this benchmark study using mass distribution data of reaction products obtained at GSI [3] in inverse kinematics (one example in FIG. 2). This step allowed us a first selection among the 10 models; in this way the first insight of the quality of the models was obtained. Then experimental mass distributions for some elements, which are interesting as RIBs will be also compared to model calculations. These data [3] have been obtained for an equivalent 0.8 or 1.0 GeV proton beam, which is approximately the proposed projectile energy. We note that in realistic thick targets the proton beam will be slowed down and some secondary particles are produced. Therefore, the residual nuclei production at lower energies is also important. For this reason, we also compared some excitation function calculations with the associated data obtained with γ-spectroscopy [START_REF] Gloris | [END_REF], [START_REF] Yu | Experimental and Theoretical Study of the Yields of Residual Product Nuclei Produced in Thin Targets Irradiated by 100-2600 MeV Protons[END_REF] to test the models in a wide projectile energy range.

We tried also to benchmark on thick target, but very few data exist. Fortunately the only ones which are available coincide with the proton beam energy and target material proposed for EURISOL, namely around 1 GeV for the energy and fissile material (ex.: U) or heavy nuclei (ex.: Pb) for the target. The first set of data, mass distributions of 5 elements, was obtained from ISOLDE [START_REF] Isolde Coll | Production yields of noble-gas isotopes from ISOLDE UCx/graphite targets[END_REF] (E p =1.0 and 1.4 GeV; target materials are ThC x and UC x ). The second set contains the specific activities of 28 radionuclides in different places along the thick target, from a dedicated experiment done at Dubna [START_REF] Pohorecki | [END_REF] (E projectile =660 MeV; target material is natural Pb). Quantitatively the MCNPX predictions are often within a factor of 2. See FIG. 3. [START_REF] Isolde Coll | Production yields of noble-gas isotopes from ISOLDE UCx/graphite targets[END_REF] with the reaction p(1.4 GeV)+UCx. On the right side, specific activities of 183 Re obtained with the reaction p(660 MeV)+Pb at Dubna [START_REF] Pohorecki | [END_REF]. The yellow rectangle indicates a too poor statistic in this region. More explanations on the figures are in the text.

FIG. 3. On the left side, mass distribution for Xe element obtained at ISOLDE

More information, details and references can be found in [START_REF] David | Benchmark calculations on residue production within the EURISOL DS project, Part I: thin targets[END_REF].

These models were also benchmarked on neutron emission, since residues in the case of the fission target option in Eurisol are produced via neutron-induced fission, neutrons coming from spallation reactions in a Hg convector (see section 4.). On this observable, with thin or thick target, all models give comparable good results. This study is developed in [START_REF] Rapp | Benchmark calculations on particle production within the EURISOL DS project[END_REF].

Taking into account the above conclusions including the benchmarks on residue and neutron production from thin and thick targets, the use of Isabel-Abla or INCL4-Abla within MCNPX 2.5.0 is recommended, and sometimes CEM2k also if one combines quality and running time.

Direct target

In order to estimate the in-target RIB yields in direct targets we had to study 320 different configurations of cylindrical targets, whose characteristics are summarized in TABLE I. The choice of the beam spot (Gaussian profile) σ = R/3 came from a previous study [START_REF] Thiollière | Optimization of 6 He production using W or Ta converter surrounded by BeO target assembly[END_REF][START_REF] Stora | Baseline parameters 05/27/2005 Eurisol-DS task #3[END_REF]. As σ values cannot be smaller than 3 mm [START_REF] Thiollière | Optimization of 6 He production using W or Ta converter surrounded by BeO target assembly[END_REF], the minimal target radius R was automatically equal to 9 mm. The other radii have been chosen to successively increase the target volume by a factor 2. Then, we have considered different target lengths. According to the ongoing EURISOL_DS, we have selected E values ranging from 0.5 to 2 GeV. Finally, as EURISOL targets are designed to stand up to 100 kW, we have fixed P to this maximal value. To get the production rates, CINDER'90 [START_REF] Wilson | A Manual for CINDER'90 Version C00D and Associated Codes and Data[END_REF], an evolution code, has been coupled to MCNPX. This code calculates the residues due to low energy neutrons, not given by MCNPX, and takes into account the decay of all isotopes. Since from our benchmarks (see section 2.), INCL4 and Isabel give similar results and CEM2k is a fast running code, we used 2 models that are INCL4/Abla and CEM2k.

Attempting to "optimize" our targets, according to elements or isotopes studied, we started by plotting mass and charge distribution for all targets and the two extreme energies. To simplify this work, radii (18mm) and masses (∼2kg) of the targets were fixed. An example is given on FIG. 4. Then bi-dimensional graphs were plotted to get optimal lengths and radii (FIG. 5). More details concerning this work can be read in [START_REF] Chabod | Optimization of in-target yields for RIB production Part I: direct targets[END_REF].

Fission target

FIG. 6 shows the target set-up as model build by MC code. The geometry model represents the last design variant, MAFF-like, able to accommodate 30kW load heat. The production targets are cylinders with internal axial holes, having a volume of 181 cm 3 (R ext = 1.75 cm, R int =0.4 cm, H=20 cm). A homogenized material of fissile-graphite compound (ρ = 1.88 g cm -3 ) with mass rate 1/20 and the actinide mass of 15 g was accounted in all calculations performed. For all configurations studied based on Uranium the moderator material was water while the reflector material was beryllium oxide. Only in case of Thorium used as fissile the iron was used in the model for both moderator and reflector materials. MCNPX2.5.0 was used in simulations.

Derivation of the fission yields was based mainly on recommended tabular yields for the fission products given in [START_REF] England | ENDF-349 Evaluation and Compilation of Fission Product Yields[END_REF]. These recommended yields include many nuclides which fission by neutrons at three energies: T for thermal energies (0.0253eV), F for fast energies (2 MeV) and H for 14 MeV for high energy. The yields were calculated using the formula:

€ Yield = T σf 0 ET ∫ E ( ) Φ E ( ) dE + F σf ET EF ∫ E ( ) Φ E ( ) dE + H σf EF 1GeV ∫ E ( ) Φ E ( ) dE (2)
where the integral limits were chosen as follow: ET = 5keV and EF = 5MeV.

In the analysis six target assemblies configurations were accounted: five cases based on Uranium compounds with 235 U percentages of: 99.99, 20, 3, 0.072 (natural Uranium) and 0.02 (depleted Uranium) as well as 232 Th.

The results presented are normalized to 1 mA intensity proton beam on the target to allow easily the desired scaling to other values. 232 Th is possible although it may reduce the magnitude of yields but yet it extends the yields range over a larger area of interests. As a matter of example FIG. 7 shows the distribution for Kr and Sn isotopes for the six studied target configurations. These elements are two of the most relevant ones for RIB production at EURISOL as indicated in [START_REF] Chabod | Optimization of in-target yields for RIB production Part I: direct targets[END_REF]. The differences in the yields between the six systems follow the fission rate trend previously discussed. All calculated isotopic distributions display an almost centred-Gauss like shape (as in the FIG. 7) that shows that the contribution of the 235 U in the mixture is dominant comparatively with 238 U whose fragment distribution is biased towards the neutron-rich wing of the distribution. 232 Th and 238 U present a slightly extend further to the 235 U fragment distribution to lighter and heavier mass. This explains the enhanced production yield in the range (72 < Z < 78) for Ni produced from 232 Th that reach values even higher than those obtained for 20% 235 U configuration.

Conclusion

After benchmarking the available spallation models in the transport code MCNPX2.5.0 for the EURISOL purpose, methodologies have been developed to estimate the in-target fission yields arising from both direct targets and converter-fission targets assembly of the EURISOL facility. This study provides quantitative estimates of the fission yields for a variety of isotopic distributions come out from different target systems. Obtained results represent helpful information for the beam intensity predictions at the future EURISOL facility and using effusion/diffusion efficiencies realistic beam intensities will be provided.

FIG. 2 .

 2 FIG. 2. Mass distribution of the reaction products from 238 U(1A.GeV)+p. GSI Data [3] are in black and model calculations are in red for the ten model combinations available in MCNPX2.5.0. Intra-Nuclear Cascades used are INCL4, Isabel and Bertini (the three first lines), and combined to various evaporation/fission models, which are Abla, Dresner/RAL and Dresner/ORNL (the three columns). The last graph (bottom left) is the stand-alone package CEM2k.

FIG. 4 .

 4 FIG. 4. Charge number distribution of the production rates for the nuclei produced insideAl 2 O 3 , Pb, SiC, Ta and UC 3 targets on the left and Pb, Ta and UC 3 targets on the right (R = 18 mm, M = π.ρ.R 2 .L = 2 kg, spallation model = INCL4/ABLA)

FIG. 5 .

 5 FIG. 5. 92 Kr (target = UC 3 , E = 0.5 GeV) on the left and 180 Hg (target=Pb,E=1GeV) production rate per incident proton [10 -6 s -1 ] as a function of target length L and radius R (model = INCL4/ABLA).

FIG 6 .

 6 FIG 6. Geometry model used in MCNPX simulations. Left: XY cross section cut through the model in the plan Z=0; Right: YZ cross section cut through the model in the plan X=0. Zones of the model are indicated in the graphs.

FIG. 7 .

 7 FIG. 7. Fission yields isotopic distributions for two selected products. Left: Kr and right Sn

TABLE II :

 II IN-TARGET PRODUCTION RATES FOR SOME ISOTOPES WITH DIRECT TARGETS USING INCL4/ABLAFinally we summarize in TABLE II optimal configurations for some isotopes and give the estimated production rate. It must be stressed that only in-target yields have been calculated,

no extraction efficiencies have been taken into account. All yields obtained here should be multiplied by the effusion/diffusion efficiencies to get realistic ion beam intensities.

TABLE III :

 III FISSION RATES FOR ALL TARGETS [FISSIONS s -1 mA -1 ] AND PERCENTAGE OF THE 235 U CONTRIBUTION.

  TABLE III shows clearly that for all accounted Uranium mixtures more than half of the total fission rates are due to fission of the 235 U component. The number of fissions (~10 15 /s) required by the baseline EURISOL project is thus met for the most favourable producing system ( 235 U case), since a 4 MW target corresponds to 4mA beam. The use of other Uranium mixtures or

We acknowledge the financial support of the European Community under the FP6 "Research Infrastructure Action -Structuring the European Research Area" EURISOL DS Project Contract no. 515768 RIDS . The EC is not liable for any use that can be made on the information contained herein.