Guilhem Mollon
email: guilhem.mollon@insa-lyon.fr

Periodic instationarities of granular flows in conical hoppers

Keywords: Granular flow, Periodic variations, Hopper flow, DEM

Granular flows through converging sections such as conical hoppers have been reported to be submitted to instationarities, which in certain circumstances can appear to be organized and periodic. In this paper, we explore this phenomenon by conducting discrete element modelling simulations of a 3D gravity-driven hopper flow and varying a large number of parameters such as hopper geometry and granular sample properties. Dedicated postprocessing techniques are developed and used to investigate the spatial and temporal patterns of these instationarities and to bring some understanding on the physics of this spontaneous phenomenon. Numerical results show that a clear structure appears for these instationarities, under the form of rapidly propagating waves relating variations in velocity magnitude and coordination number. While very faint, periodic variations of the sample density are also detected. The parametric study reveals that the self-organization of these variations requires a narrow set of conditions in terms of hopper geometry and intergranular contact friction coefficient.

Introduction

Hopper flow is a common problem in physics of granular flows [START_REF] Forterre | Flows of dense granular media[END_REF]. It consists in the gravity-driven flow of a granular assembly through a converging channel, and the classic hourglass problem is a good example of this class of situations. Besides its theoretical interest, this problem also has a large number of practical applications since converging flows are used in a wide variety of industries manipulating powders and grains [START_REF] Cleary | DEM modelling of industrial granular flows: 3D case studies and the effect of particle shape on discharge[END_REF][START_REF] Wan | Influence of geometrical and material parameters on flow rate in simplified ADS dense granular-flow target: a preliminary study[END_REF][START_REF] Jia | Fluctuation and arching formation of very dense and slow pebble flow in a silo bed[END_REF][START_REF] Pizette | DEM GPU studies of industrial scle particle simulations for granular flow civil engineering applications[END_REF]. It has therefore attracted a large attention from the scientific community, both from the practical and academic viewpoints. Theoretical [START_REF] Nguyen | Gravity flow of granular materials in conical hoppers[END_REF][START_REF] Bazant | A theory of cooperative diffusion in dense granular flows[END_REF][START_REF] Hendy | Instabilities in granular flows[END_REF][START_REF] Sun | Radial hopper flow prediction using constitutive model with microstructure evolution[END_REF], experimental [START_REF] Michalowski | Flow of granular material through a plane hopper[END_REF][START_REF] Baxter | Pattern formation in flowing sand[END_REF][START_REF] Choi | Velocity profile of granular flows inside silos and hoppers[END_REF][START_REF] Gardel | Forcevelocity correlations in a dense, collisional, granular flow[END_REF][START_REF] Gardel | Dynamical fluctuations in dense granular flows[END_REF][START_REF] Gentzler | Measurement of velocity and density profiles in discharging conical hoppers by NMR imaging[END_REF][START_REF] Vivanco | Dynamical arching in a two dimensional granular flow[END_REF][START_REF] Villagran Olivares | Towards a one parameter equation for a silo discharging model with inclined outlets[END_REF] and numerical [START_REF] Ristow | Density patterns in two-dimensional hoppers[END_REF][START_REF] Potapov | Computer simulation of hopper flow[END_REF][START_REF] Mollon | Characterization of fluctuations in granular hopper flow[END_REF][START_REF] Magalhães | Analysis of the velocity field of granular hopper flow[END_REF][START_REF] Mollon | Granular flow fluctuations in a conical hopper[END_REF] studies have been reported using a number of different tools.

The most common quantity of interest that is extracted from such studies is the hopper discharge rate, i.e. the amount of matter passing through the hopper outlet in a certain amount of time. Contemporary studies use the hopper flow configuration to investigate such complex topics as avalanching [START_REF] Wang | Velocity profiles of avalanches during hopper discharge[END_REF], self-structuration of the granular fabric triggered by the converging flow [START_REF] Guo | Enhancing the linear flow of fine granules through the addition of elongated particles[END_REF], shape-induced flow regimes [START_REF] Szabo | Flow of anisometric particles in a quasitwo-dimensional hopper[END_REF], or phase change (i.e. jamming) of the granular sample [START_REF] Tang | How granular materials jam in a hopper[END_REF]. Despite the simplicity of its geometry, the hopper flow gives rise to fairly complex behaviors, and one of the most mysterious is the spontaneous emergence of periodic flow instationarities. Such instationarities, sometimes called fluctuations, have been reported experimentally [START_REF] Baxter | Pattern formation in flowing sand[END_REF][START_REF] Gardel | Forcevelocity correlations in a dense, collisional, granular flow[END_REF][START_REF] Gardel | Dynamical fluctuations in dense granular flows[END_REF][START_REF] Vivanco | Dynamical arching in a two dimensional granular flow[END_REF] and numerically [START_REF] Ristow | Density patterns in two-dimensional hoppers[END_REF][START_REF] Mollon | Characterization of fluctuations in granular hopper flow[END_REF][START_REF] Mollon | Granular flow fluctuations in a conical hopper[END_REF].

In particular, a numerical study accounting for realistic sand grains shapes but limited to 2D kinematics [START_REF] Mollon | Characterization of fluctuations in granular hopper flow[END_REF] showed that these periodic instationarities could be described as waves propagating upstream, triggered by the build-up of force chains from the hopper walls. In the present study, we enrich this numerical framework by adding a third dimension and deploying appropriate statistical methods to extract knowledge on this phenomenon. Section 2 presents the numerical framework, including the postprocessing techniques, and describes the main flow parameters for a reference case. Section 3 provides a detailed analysis of the instationarities detected for this reference case, and Sect. 4 adds a parametric study involving a large number of simulations and systematically investigating several parameters of the system. Discussion of the results and proposals for future research avenues are presented in Sect. 5.

Numerical framework

Simulations principles and parameters

The Discrete Element Method (DEM, [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF]) has now become a very common investigation tool in granular science, both for academics and practitioners, and is used in a large variety of scientific fields such as geophysics [START_REF] Mollon | Discrete modelling of rock avalanches: sensitivity to block and slope geometries[END_REF] or tribology [START_REF] Mollon | A numerical framework for discrete modelling of friction and wear using Voronoi polyhedrons[END_REF]. It is based on an explicit integration in time of the Newtonian dynamics of each one of the bodies composing the granular assembly. In addition to a possible gravity field, these bodies are submitted to interaction forces that arise when they enter into contact with each other. Various phenomenological contact laws can be used, but we will stick in the present study to the standard contact model of DEM. It includes contact stiffnesses in the normal (k n ) and tangential ( k t ) directions, which can be seen either as purely numerical regularization parameters or as a proxy for the local deformations of the contacting bodies. In the latter case, non-linear stiffnesses reproducing Hetrz-Mindlin theory can be used. In the present study, however we remain in the first interpretation of these parameters, and use k n = k t = constant . The standard model also includes two modes of energy dissipation: a classical dashpot introducing a damping in the normal direction, and a Coulomb friction coefficient in the tangential direction.

In this paper we perform DEM simulations of perfectly conical hopper discharge, with various sample properties and hopper geometries. All the granular samples are composed of 50,000 spherical grains with a mean diameter D = 1 . Since the study is dimensionless, this quantity will be used as a length unit. The size distribution of the grains around this mean value is homogeneous, with coefficients of variation COV D (i.e. ratio of the standard deviation to the mean value) of D varying from 0% (perfectly monodis- perse sample) to 80% (strongly polydisperse sample). The opening diameter L of the hopper varies between 6D and 14D , and the hopper angle to the horizontal varies from 40 • to 80 • (Fig. 1). The gravity is kept constant and equal to 10, and the density of the grains is equal to 1. Hence, a given simulation is defined by 6 parameters, and the reference simulation corresponds to COV D = 40% , L = 10D , = 70 • , k n = 10 5 , = 0.4 , and = 0.2 . The results of this simulation are detailed in the next section. The time step is constant and equal to 10 -5 time units for all the simulations reported here. The only characteristic times of the system are related to free fall and to single-grain harmonic oscillations, and are respectively equal to:

T g = √ g∕D ≈ 0.316 T k = √ 6k n ∕ D 3 ≈ 0.00229

Postprocessing techniques

The main quantity of interest when studying hopper flow is the discharge rate, since it is of primary importance regarding the practical use of such a system for engineering applications. Figure 2 shows the evolution in time of the proportion of granular material passed through the hopper outlet, as predicted by the reference simulation. The passing rate is constant, in good agreement with existing theories and observations. This is also the case for all the other simulations described in this work. Since this paper is primarily interested in granular flow instationarities, a dedicated postprocessing technique is needed in order to follow the evolution of quantities which vary both in space and time. First, in order to avoid a drift in the average values of local flow parameters, the focus is put on a short time interval used as an observation window. As shown in Fig. 2, this interval runs from 5 to 7 time units after opening of the hopper, in order to make sure that the discharge flow has reached a steady state (in some simulations described in the next sections, a different starting time was chosen for the same reason). This time interval is short enough to consider that the mean flow parameters do not change much, but long enough to observe a statistically significant number of fluctuations.

A major difficulty is the fact that relevant quantities are attached either to one grain (velocity and coordination number, i.e. number of contacting neighbors) or to a collection of grains (solid fraction). This is a logical consequence of the Lagrangian framework used to describe the problem, but renders direct interpretations difficult because these grains (or groups of grains) are in constant motion. This issue was already tackled in several papers with various mathematical approaches [START_REF] Goldhirsch | Stress, stress asymmetry and couple stress: from discrete particles to continuous fields[END_REF][START_REF] Weinhart | Closure relations for shallow granular flows from particle simulations[END_REF]. In this paper, a specific technique is used in order to transfer these quantities to fixed coordinates and study the flow from an Eulerian perspective. Grain-attached quantities (velocity and coordination number) are interpolated at fixed locations between grains using moving least squares [START_REF] Fasshauer | Meshfree approximation methods with matlab[END_REF], in a spherical radius equal to 5 distance units around the desired location (it is important to notice that a certain amount of spatial averaging is included in the process). And a local solid fraction at a given fixed location (x, y, z) of space is computed by a weighted averag- ing technique which complies with the principle of mass conservation:

In this expression, V i corresponds to the volume of the grain i , r i (x, y, z) corresponds to the distance from the center of this grain to the location (x, y, z) , and the summation is executed on all the grains i located in a sphere of radius d = 8D around this location, which thus controls the cut-off interpolation distance. These post-processing techniques are used to extract and analyze the evolutions of velocity, coordination number and solid fraction at various locations, with a time step of 0.002 time units. A 3D Eulerian grid is used for that purpose, as shown in Fig. 2: 330 points are defined in a given vertical half-plane containing the symmetry axis of the hopper, and 8 such half-planes are used in order to cover the whole system.

F s (x, y, z) = d 2 3∕2
⋅ ∑ V i e -r i (x,y,z)∕d

Mean fields

Average fields of vertical and radial velocities, of coordination number and of solid fraction are plotted in Fig. 3 (in one half-plane, after averaging in the circumferential direction), as well as the corresponding time-related standard deviations of these quantities, for the reference simulation. In terms of time-average, the radial component of the velocity field is small compared to the vertical one, which progressively increases in norm when getting closer to the hopper outlet.

The coordination number and the solid fraction are almost constant in the hopper (close to 3.5 and 0.45 respectively), except in the neighborhood of the inlet and of the outlet of the hopper where these values drop because of reduced local confining stress. These behaviors are consistent with common hopper flow observations (e.g. [START_REF] Vivanco | Dynamical arching in a two dimensional granular flow[END_REF][START_REF] Mollon | Characterization of fluctuations in granular hopper flow[END_REF]).

The time-variability of the velocity field presents some interesting spatial patterns. Two areas are observed where velocity has important variations in time: the area located close to the outlet (because of local agitation related to the change of flow regime when the grains are released from the hopper confinement) and a broad area located in the upper half of the hopper, starting about 25 -30D above the outlet and extending up to the inlet. There is also a large time variability of the coordination number in the same area, while the solid fraction variability is very low in the whole sample. 3 Analysis of the instationarities

Detection

In order to describe quantitatively the flow instationarities, we perform the same postprocessing approach on a vertical profile corresponding to the axis of revolution of the hopper, in the reference case. The sampling time-step remains equal to 0.002 time units, and the sampling step in space along the vertical axis is equal to 0.01D . We therefore obtain data under the form Q(t, Y) where Q is velocity magnitude, coordination number, or solid fraction, and where t and Y denote time (with a zero corresponding by convention to the beginning of the observation time window, for the remainder of the present study) and vertical position (with a zero corresponding to the hopper outlet).

These data are represented as space-time diagrams in the left-hand column of Fig. 4. Some organized fluctuations clearly appear in the diagrams of velocity and coordination number (but not in the diagram of solid fraction). It is however difficult to read and quantify them, because they are superimposed to a mean field which is spatially varying (as well as its time-variability, as shown in Fig. 3). This is especially the case for the velocity, which progressively increases when getting closer to the outlet. For this reason, quantities are locally normalized, using the following formula:

In this expression, Q(t, Y) is the normalized quantity, Q(t, Y) is the initial quantity, Q (Y) is the time-average of Q at the location Y , and Q (Y) is the corresponding time-related standard deviation. Normalized velocity, coordination number, and solid fraction are plotted on space-time diagrams on the right-hand column of Fig. 4. Fluctuation patterns are much clearer in that normalized form: for both velocity and coordination number, they take the form of periodic oscillations around the mean value, which seem to be strongly synchronized in the vertical directions. More precisely, the lower part of the hopper does not seem to be subjected to any organized pattern, but above a certain height (which is close to Y = 15D for the velocity and Y = 10D for the coordination number), fluctuations acquire a structure which remains the same up to the hopper inlet. The case of solid fraction is more complex: normalized solid fraction seems also to be submitted to periodic and synchronized oscillations in the upper part of the hopper, but these oscillations are fainter and are blurred by noisy descending patterns which may be related to the downward motion of undisturbed grains structures within the granular flow.

Q(t, Y) = Q(t, Y) -𝜇 Q (Y) 𝜎 Q (Y)
It is interesting to notice that the periods of oscillation of the three considered quantities are the same, indicating that these oscillations are part of the same physical phenomenon. In similar studies in 2D [START_REF] Vivanco | Dynamical arching in a two dimensional granular flow[END_REF][START_REF] Mollon | Characterization of fluctuations in granular hopper flow[END_REF], it was observed that such oscillations could be interpreted as waves emanating from the hopper outlet and travelling upwards in the hopper. The velocity of these waves was reported to be much greater than that of the flowing grains. This interpretation is debatable here since the oscillations described in Fig. 4 seem to be synchronized in space. At a first glance, indeed, they look like stationary oscillations and do not seem to exhibit any propagation pattern. However, this could be related to such waves travelling too fast (with respect to their time-period) for this propagation to be detectable visually, because the hopper is not high enough to allow for a sufficient propagation distance. A preliminary analysis of similar 3D simulations [START_REF] Mollon | Granular flow fluctuations in a conical hopper[END_REF] indeed revealed that the use of a lower value of the contact stiffness between grains was able to reduce the period of the oscillations and to slow-down their propagation, allowing a clear observation of the wave upwards travel. We can thus hypothesize that the oscillations which are described here result from the same phenomenon than in 2D [START_REF] Mollon | Characterization of fluctuations in granular hopper flow[END_REF], but with waves travelling much faster.

Patterns in time

Figure 5 shows an illustrative snapshot of the spatial fields of normalized velocity and coordination number in a vertical plane of symmetry of the hopper. This visualization reveals that such instantaneous fields exhibit some clear trends (for example in the vertical direction) but an important spatial variability. Hence, the extraction of meaningful patterns requires some specific procedures.

For this purpose, we explore the autocorrelation structure of the time variations of the normalized quantities in Fig. 6. At a given location Y on the revolution axis of the hopper, a given quantity Q(t, Y) varies in time, and its autocorrela- tion function A Q (Δt, Y) can be computed as a function of the phase shift Δt by computing the correlation between Q(t, Y) and Q(t +Δt, Y) . Maps of this autocorrelation are plotted in Fig. 6 for velocity, coordination number, and solid fraction and reveal interesting patterns. For both velocity and coordination number, a very clear autocorrelation peak occurs for a phase shift Δt = T = 0.192 time units , indicating a strong periodicity. This is only the case, however, beyond a certain height above the hopper outlet (close to Y = 15D for the velocity and Y = 10D for the coordination number), in good accordance with the observations of Fig. 4. Below this height, the correlation structure disappears. Besides, a similar map of solid fraction does not show any clear autocorrelation peak, indicating that the faint periodic signal observed in Fig. 4 is not strong enough to be revealed by this procedure. It is interesting to notice that the measured In order to describe more precisely the waveform of the detected oscillations, we plot in Fig. 7 the evolution in time of the normalized velocity magnitude, coordination number, and solid fraction at selected illustrative locations. More specifically, we focus on eight different locations in the hopper which share the same Y-coordinate and the same distance to the hopper axis, but in eight different radial directions (corresponding to the eight half-planes described in Fig. 2). As expected, these signals are very well synchronized (over about ten periods) for velocity and coordination number but do not show any organized pattern for the solid fraction. These signals are then stacked (i.e. averaged on the 8 locations and on 10 successive periods, using the main period value T = 0.192 obtained in the previous paragraph), as is often done in seismology to de-noise earthquakes waveforms [START_REF] Grigoli | Automated microseismic event location using master-event waveform stacking[END_REF]. The purpose of this operation is to focus on the systematic variations of the quantities that follow the main oscillation period while getting rid (by averaging) of all the other non-systematic variations. The resulting waveforms are plotted in the lower-part of Fig. 7 for the three quantities (velocity, coordination number, and solid fraction). Two identical periods of this waveform are shown to ease the reading.

Several interesting observations can be extracted from these curves. It appears that the velocity oscillations are quasi-harmonic, with an evolution which is very symmetric around the average value (equal to 0 for such normalized quantities). The coordination number waveform, on the other hand, is strongly skewed and far from being harmonic. The main pattern that this quantity follows in time is divided in a rather long period of values moderately larger than the average (up to + 1 standard deviation), followed by a shorter but stronger drop below the average (up to -2 standard deviations). The stacking operation also reveals a systematic evolution of the solid fraction which could not be inferred from autocorrelation patterns (Fig. 6) or from the raw timesignals (upper part of Fig. 7). While noisier than the other quantities, the waveform of the solid fraction is very clear, and consists in a symmetrical signal, somewhere between harmonic and triangular waveforms.

Solid fraction and coordination number waveforms are in phase, indicating that these two quantities are directly related. As could be expected, a larger solid fraction correlates with an increase in the number of contacts, and vice versa. Their waveforms are however different because of the intrinsically asymmetrical character of mechanical contacts. Indeed, since the stiffness of the contacts is pretty large, the interpenetration distance of contacting grains is always quite small. A minor decrease of the solid fraction is hence sufficient to increase a little bit the distance between the centers of the grains, which leads to the loss of a large number of contacts. Conversely, when the solid fraction increases, some contacts can be created if the initial gap separating the grains is small enough for them to eventually touch, but this situation is made rarer by the large stiffness of the contact model. This is the reason why the coordination number oscillations are skewed towards the low values while the solid fraction follows a more harmonic evolution.

In contrast, the waveforms of velocity and coordination number are shifted in time: the velocity oscillations are delayed by about 0.3T with respect to the variations of coor- dination number. Broadly speaking, it appears that velocity magnitude decreases when the coordination number is above its average, and increases when it is below its average. Likewise, coordination number decreases when velocity is below its average, and increases when velocity is above its average.

Patterns in space

This stacking technique is further applied to all locations in the hopper half-plane in order to analyze the spatial patterns of evolution of the three quantities of interest, and the resulting fields are plotted for one period T in Fig. 8.

The lower part of the hopper appears to have no significant variation in this representation. It does not mean that the fields do not vary in time at these locations, it means that staking the signals on a T-period averages them and thus that the signals do not have any T-periodicity in this area. The upper row of Fig. 8 shows that the velocity oscillations only appear above a height of about 15 -20D (con- firming the observations of Figs. 6 and7), and that they are very homogeneous. This area can be described as the organized fluctuations area. The same observation applies to the coordination number, although it clearly appears that organized fluctuations do start at a lower height of 10 -15D , and that they are very slightly less homogeneous This view is confirmed by the lower part of Fig. 9, which shows maps of wave delay, provided by the stacked signals at every locations of the half-plane. More specifically, it shows a map of the arrival time of the maximum velocity (lower-left part) and of the minimum velocity (lower-right part). The time delay between two successive lines in the maps is constant and equal to 0.0005 time units. It clearly appears that the velocity increase initiates Fig. 7 Pattern extraction for the velocity, the coordination number, and the solid fraction, at an illustrative location in the hopper. Upper-part: time-signals of normalized quantities obtained on eight similar locations of the hopper (same height and distance to the axis, different radial directions as shown in Fig. 2); lower-part: stacked signals for the three quantities on the lateral wall of the hopper, propagates inwards, and then upwards. The velocity decrease, on the contrary, initiates at the walls but also at the inlet and below the organized oscillation area, and then propagates inwards. It also propagates about three times faster. These observations confirm that the detected instationarities are indeed waves propagating in the hopper, and not only stationary vibrations. In the considered geometry, the propagation time remains however very small when compared to the main period of the oscillations.

Parametric study

Influence of the particle size distribution

In this section, we investigate the influence of the size variability of the granular sample on the flow instationarities. This variability is quantified by the coefficient of variation of the particles diameters COV D , i.e. the ratio of the stand- ard deviation to the mean value of the diameter (which is always equal to 1 length unit in this study). COV D takes the Fig. 8 A typical oscillation cycle of period T starting at time t 0 (Fig. 7), obtained by stacking for the three quantities of interest in a hopper halfplane values 0%, 20%, 40% (reference case), 60%, and 80%, while all the other simulation parameters are kept constant and equal to those of the reference case. Results are presented in Fig. 10, and include (1) the discharge rate of the hopper flow during the whole discharge (this rate is found to be constant throughout the whole simulation, for all simulations presented in this work), ( 2) the average and standard deviation of the velocity magnitude at a height Y = 40D (i.e. in the organized oscillation area) and the corresponding COV, (3) the same data for coordination number, and (4) the maximum velocity autocorrelation at Y = 40D and the corresponding time-period (obtained in a manner similar to Fig. 6).

These results reveal that the flow rate through the hopper outlet is unaffected by the size distribution of the grains, as well as the average velocity in the organized oscillation area, while the average coordination number slightly increases with COV D . The intensity of the instationarities (i.e. the coefficient of variation of these quantities) seems to decrease when increasing the dispersion of the grain sizes in the sample, but the coherence of the periodicity (quantified by the value of the maximum autocorrelation number) and their time-period do not seem to be modified by this parameter. Hence, polydispersity seems to limit the instationarities, but neither to change their nature nor to prevent them. Following this trend, we can however speculate that a very strongly polydisperse sample (e.g. with a fractal distribution) may only experience limited amount of flow oscillations, but it would have to be confirmed by further studies.

Influence of the hopper geometry

Figure 11 presents the influence of the diameter of the hopper outlet, with values varying from 6D to 14D (all the other parameters being constant and equal to those of the reference experiment). Flow rate and average velocity both increase quite strongly with the outlet diameter, as predicted by Beverloo law, and the average coordination number is rather small (3.3-3.4) for intermediate values of the outlet diameter ( 8 -10D ) but larger (3.8-4.2) for extreme values.

We find the same kind of trend when analyzing the intensity of the instationarities (COV of velocity and Fig. 9 Typical velocity signals (obtained by stacking) at several locations in a hopper half-plane (with a small vertical shift to ease the reading), maps of signal delays, and corresponding propagation directions; upper-left: stacked signals on a vertical profile; upper-right: stacked signals on a radial profile; lower-left: arrival time of the maximum of velocity; lower-right: arrival time of the minimum velocity coordination number), the quality of their periodicity structure (maximum autocorrelation value) and their timeperiod. It is clear that hopper outlet diameters of 8 -10D favor the appearance and development of the organized flow fluctuations. For these simulations, the instationarities are more intense (with coefficients of variation of 30-50% for velocity and 30-40% for coordination number), and with clearer patterns (maximum autocorrelation close to 1), with a constant period. In contrast, when considering very small ( 6D ) or large ( ≥ 12D ) outlet diam- eters, the oscillations get less intense and with a much less clear periodicity. This is especially the case for large diameters, with coefficients of variation smaller than 10% for the velocity and the coordination number, a maximum autocorrelation lower than 0.8 (and as low as 0.5 in the case of 14D ), and a decreasing time-period. It is likely that considering even larger outlet diameters would lead these organized fluctuations to vanish.

The influence of the hopper opening angle is presented in Fig. 12, from the case of a very flat hopper (angle to the horizontal = 40 • ) to a very steep one ( = 80 • ). Since the general geometry of the hopper changes in each simulation, local quantities (i.e. mean and standard deviation of velocity and coordination number) are computed at a varying height located 10D below the free surface at the hopper inlet, which corresponds to the location where organized fluctuations are the most prone to be observed.

Influence of the contact model

In a preliminary study performed in the same numerical framework [START_REF] Mollon | Granular flow fluctuations in a conical hopper[END_REF], the effect of the contact stiffness k n = k t was investigated. It was found that this parameter did not have any influence on the discharge rate, but that the periodicity of the flow oscillations was proportional to the inverse square root of this stiffness. It clearly indicated that the development of such oscillations relied on the elasticity of the contacts between the grains.

In this paragraph, we evaluate the influence of the two other parameters of the standard contact model of DEM: friction and damping coefficients. Figure 13 presents numerical results obtained for friction coefficients varying from = 0.0 to = 0.8 . The hopper discharge rate slightly decreases with an increasing friction, as well as the average velocity and coordination number. There is however a clear change of oscillation regime between low and high friction values. For friction coefficients lower than 0.4, instationarities are very faint (coefficients of variation close to 5% for velocity and coordination number) and only poorly periodic. In the case of a null friction, the periodic structure completely vanishes. For larger friction coefficients, however, the intensity and coherence of the oscillations are high and rather constant. It is interesting to notice that large friction coefficients lead to larger time periods of the oscillations.

The influence of the damping coefficient is detailed in Fig. 14, for damping values ranging from = 0.0 to = 0.4 (a value of 1.0 corresponding to critical damping). It appears that the influence of this parameter on the system is rather limited, except for a decrease of the intensity of the coordination number oscillations with a larger damping. For all the

Discussion and conclusion

The observations gathered during this numerical campaign reveal that 3D granular hopper flows are very prone to organized instationarities. They appear to occur as intense variations of the velocity field and of the contact structure of the granular sample, which are chaotic in the neighborhood of the hopper outlet but gain a clear periodicity and structure above a certain height. They take the shape of rapidly propagating waves originating from the hopper outlet and lateral walls and propagating upwards and inwards, leading to almost synchronized oscillations of the flow parameters in the upper part of the hopper. While velocity variations are quasi-harmonic, variations in the coordination number are more skewed, with short events of important contacts loss followed by longer stages of reconstruction of the contact network. The granular density of the sample, represented here by the solid fraction, is only submitted to very limited variations in amplitude, but can nevertheless be detected by appropriate procedures and appear to be perfectly in phase with the coordination number signal.

The time delay between the velocity and coordination number variations and the propagation patterns of the waves transporting them draw a general pictures of the instationarities that follow several successive stages:

1. Force chains build up from the hopper outlet and from its lateral walls, because of the converging character of the flow. This leads to an increase in the coordination num- The systematic parametric study described in Sect. 4 suggests that some specific conditions are necessary for the initiation and continuance of such organized and periodic instationarities. A specific range of hopper outlet diameter (between 8D and 10D ) and of hopper open- ing angle (between 60° and 70° with the horizontal) are needed, otherwise the fluctuations do not gain intensity and periodic structure. A minimum value of the coefficient of friction of 0.4 is also necessary. The polydispersity of The reasons for these restricted conditions are still unclear, in particular because there is no apparent correlation between the discharge rate and the properties of the instationarities. For example, the change of oscillation regime is very clear between the cases = 0.2 and = 0.4 (Fig. 13), but no such drastic change is observed in terms of the discharge rate, which follows the same trend as for larger values of . Hence, a simplistic explanation relating the appearance of the fluctuations to the larger difficulty of the grains to flow through the hopper (which would lead them to collide more and would induce vibrations) does not seem to hold. Same observations can be made for the hopper outlet diameter and opening angle. This is confirmed by Fig. 15, which shows vertical profiles of the coefficient of variation of the velocity magnitude and of the coordination number along the revolution axis of the hopper, for the reference simulation, for the case of hopper opening diameter of 6D , and for the case of a fric- tion angle of = 0.2 . It clearly appears that the behaviors of these three systems in terms of flow variability are very similar in the lower part of the hopper, i.e. below a height of 10 -15D . It is only above this height that the reference sim- ulation starts to exhibit more intense instationarities, while the other two simulations do not. It would thus appear that organized fluctuations are not related to phenomena occurring in the close neighborhood of the outlet, such as the collisional regime transition. It can be speculated that their ability to self-organize may result from a complex interplay between the hopper geometry and the ability of the granular material (enhanced by a higher contact friction) to form strong and persistent force chains and thus to coordinate grains motions on a longer range. This interplay is yet to be understood. Further studies will focus on this assumption, and will explore the implications of a more realistic (i.e. non-linear) contact stiffness, of larger polydispersity, and of more complex grains shapes [START_REF] Mollon | 3D generation of realistic granular samples based on random fields theory and Fourier shape descriptors[END_REF].

Fig. 1

 1 Fig. 1 Snapshot of particle velocities during hopper discharge

Fig. 2 Fig. 3

 23 Fig. 2 Left: Proportion of passed material in time in the reference simulation, and timewindow used for the observation of flow instationarities; right: 3D Eulerian grid used to map flow parameters in time

Fig. 4

 4 Fig. 4 Space-time diagram of the velocity magnitude, of the coordination number, and of the solid fraction, on the vertical revolution axis of the hopper; left: raw quantities; right: normalized quantities

Fig. 5

 5 Fig. 5 Snapshot of normalized flow parameters in a symmetry plane of the hopper

Fig. 6

 6 Fig. 6 Autocorrelation in time of normalized flow quantities on the revolution axis of the hopper, presented as functions of location Y and of phase shift; upper-left: Velocity autocorrelation map; upper-

  Figure 12 reveals that steep hoppers tend to increase the flow rate, but that fluctuations are favored by intermediate values of the opening angle. Flat hoppers ( = 40 -50 • ) indeed lead to very limited insta- tionarities (COV smaller that 8% for velocity and coordination number) with a much less clear periodic structure. The case = 40 • even leads to a complete disappearance of any kind of periodic fluctuations of the flow parameters. Similar observations are made for the steepest case ( = 80 • ). It is nevertheless interesting to observe that the shape of the

Fig. 10

 10 Fig. 10 Influence of the coefficient of variation of D ; upper-left: discharge rate; upper-right: velocity magnitude at Y = 40D (mean value ± one SD), and corresponding coefficient of variation; lower-

Fig. 11

 11 Fig. 11 Influence of the outlet diameter; upper-left: discharge rate; upper-right: velocity magnitude at Y = 40D (mean value ± one SD), and corresponding coefficient of variation; lower-left: same plot for

Fig. 12

 12 Fig. 12 Influence of the hopper angle; upper-left: discharge rate; upper-right: velocity magnitude at a distance of 10 below the hopper inlet (mean value ± one SD), and corresponding coefficient of varia-

Fig. 13

 13 Fig.[START_REF] Gardel | Forcevelocity correlations in a dense, collisional, granular flow[END_REF] Influence of the contact friction coefficient; upper-left: discharge rate; upper-right: velocity magnitude at Y = 40D (mean value ± one SD), and corresponding coefficient of variation; lower-

Fig. 14

 14 Fig. 14 Influence of the contact damping coefficient; upper-left: discharge rate; upper-right: velocity magnitude at Y = 40D (mean value ± one SD), and corresponding coefficient of variation; lower-

Fig. 15

 15 Fig.[START_REF] Gentzler | Measurement of velocity and density profiles in discharging conical hoppers by NMR imaging[END_REF] Standard deviations of velocity and coordination number as functions of the location Y on the revolution axis of the hopper, for the reference simulation, the case of an opening diameter of 6D , and the case of a friction coefficient = 0.2


Olivier Bouillanne

Guilhem Mollon
email: guilhem.mollon@insa-lyon.fr

Aurélien Saulot

Sylvie Descartes

Nathalie Serres

Guillaume Chassaing

Karim Demmou

How vorticity and agglomeration control shear strength in soft cohesive granular flows

Keywords: Soft grains, Sheared flow, Agglomeration, Vorticity, Tribological contact, Third body

Deformable granular flows present complex kinematics. These materials can have various flow regimes: plastic, agglomerated, rigid-like granular flow, etc. In this paper, a multibody meshfree model is used to investigate the consequences of cohesion, stiffness, and viscosity of the particles on their collective sheared flows in tribological contacts. An approach derived from fluid mechanics postprocessing tools, based on vortex detection, is employed to understand the links between these parameters and the emerging friction coefficient of the sheared interface. These results pave the way to complete kinematic studies of third body simulations in tribological contacts.

Introduction

Tribology and fretting

Tribology is the science that describes the interaction of the surfaces of two solid bodies in contact and relative motion. This contact generates friction and can lead to various mechanical phenomena, such as wear and fatigue. Friction is at the center of many practical topics, from moving cut stones by the ancient Egyptians [START_REF] Dowson | History of Tribology[END_REF] to brushing teeth [START_REF] Descartes | Tribological study of oral care silica[END_REF]. Mechanical systems can be lubricated or dry. While the understanding of lubricated contacts phenomenology is now fairly accurate, that of dry contacts remains fragmentary.

An approach popularized by Godet [START_REF] Godet | The third-body approach: a mechanical view of wear[END_REF] and then Berthier [START_REF] Berthier | Velocity accommodation in fretting[END_REF][START_REF] Berthier | Fretting fatigue and fretting wear[END_REF] is to consider that, at the interface between two rubbing bodies, there is a third body. This third body, composed of a combination of wear debris from the first bodies and particles from the outside, is located between the surfaces and allows loads to be transmitted while accommodating the relative velocity of the two rubbing bodies. This view is supported by countless experimental characterizations of worn surfaces. Hereinafter, these two rubbing bodies will be called "first bodies", following the conventional tribological terminology.

Fretting is a tribological phenomenon in which two surfaces in contact undergo oscillatory movements of amplitude smaller than the size of the contact. This fretting can generate fatigue or wear, or a mixture of both. Many mechanical parts undergo fretting: splines, bolted or riveted assemblies [START_REF] Kounoudji | Role of third body on bolted joints' self-loosening[END_REF], etc. This affects many fields of engineering: aeronautics, nuclear, biomedical, civil engineering, etc. The concept of third body has led to a better understanding of fretting, especially following advances made by the fretting maps [START_REF] Vingsbo | On fretting maps[END_REF] and the study of the different fretting regimes (stick condition, stick-slip condition, and gross-slip condition) [START_REF] Vincent | Mechanics and materials in fretting[END_REF].

The third body may protect or damage surfaces and can show various aspects, as can be seen in Fig. 1. These SEM views are obtained as a result of two titanium alloy (Ti-6Al-4 V) surfaces rubbing against each other in a frettingtest tribometer, i.e. following a small amplitude reciprocating sliding motion during a large number of cycles. As in any such test, a third body is visible. On the left-hand micrograph, various third body behaviors are detectable: granular on the left, and rather plastic in the center. On the right-hand micrograph, the third body is granular and large agglomerated wear debris are visible.

Experimental approaches can reveal this third body after the test, but traditional mechanical models tend to neglect the impact that the third body can have on wear. The main reason for this is the fact that state-of-the-art measurement techniques are not able to resolve in time and space the variability of the local third body behaviors and their influence on local friction and wear. There is therefore a lack of knowledge of such local effects, and the best way to circumvent this experimental bottleneck is to employ numerical simulations.

Conventional numerical methods of modeling wear are descriptive and hardly predictive. This is due to the multitude of factors influencing interface behavior: the amplitude of the movement between the two surfaces, the pressure distribution on these surfaces, the number of accumulated cycles and their frequency, the properties of the materials and the condition of their surfaces, the contact temperature, the atmospheric environment… Collins [START_REF] Collins | Failure of Materials in Mechanical Design[END_REF] thus counts about fifty variables that play a role in the fretting process. The relationship between the third body and wear is established but not yet totally predictable [START_REF] Arnaud | Modeling adhesive and abrasive wear phenomena in fretting interfaces: a multiphysics approach coupling friction energy, third body and contact oxygenation concepts[END_REF][START_REF] Baydoun | Modeling contact size effect on fretting wear: a combined contact oxygenation-third body approach[END_REF].

Surface degradation leads to wear after a while, especially if the particles that have become detached from the surface come out of contact permanently. If these particles remain within the contact, they form the third body and can help to limit this wear. The stresses applied by the flow of the third body in the interface are probably one of the causes of surface degradation. It is therefore necessary to characterize in detail the stresses applied to the first bodies by the third body. Such forces are obviously related to the flow regime of the third body in the interface, especially on its shear strength and on its spatial and temporal variations. Kinematic indicators will allow us to characterize these regimes and to understand the link between the third body and the stresses, then with the degradation, and finally with the wear.

The kinematics of third body flows must therefore be better understood, to grasp the interaction between the third body and stresses. Numerical [START_REF] Wang | Strain dependent vorticity in sheared granular media[END_REF] and experimental [START_REF] Forterre | Stability analysis of rapid granular chute flows: formation of longitudinal vortices[END_REF][START_REF] Abedi | Vortex formation and dissolution in sheared sands[END_REF][START_REF] Rognon | A circulation-based method for detecting vortices in granular materials[END_REF] approaches to calculate the vorticity of a granular flow have already been proposed, but they were limited to rigid granular materials. The terms agglomerate and vortex are complementary. The term agglomerate refers to a set of particles moving in a coordinated manner, while a vortex refers to a set of particles rotating around a certain geometrical point. Hereafter, in this study, the term agglomerates or structures will be used, because the vortices studied here have much longer duration than those described in the existing granular physics literature.

Numerical tribology

Different numerical approaches are used to describe the tribological contacts. The first one is based on finite elements. To model the wear track, one way to do this is to calculate the stresses undergone by the surfaces and to deduce a displacement of the nodes composing the mesh of these solids. [START_REF] Ding | The effect of slip regime on fretting wear-induced stress evolution[END_REF]. Another way, more radical, consists of recalculating a complete mesh as the wear progresses [START_REF] Basseville | An evaluation of the competition between wear and crack initiation in fretting conditions for Ti-6Al-4V alloy[END_REF]. The "wearbox" numerical tool, on the other hand, relies on a measurement of the amount of energy dissipated in contact to apply a deformation to the mesh [START_REF] Mary | Numerical prediction of fretting contact durability using energy wear approach: optimisation of finiteelement model[END_REF][START_REF] Paulin | Finite element modelling of fretting wear surface evolution: application to a Ti-6A1-4V contact[END_REF]. The shortcoming of these methods is that they do not consider the debris bed present in the interface, i.e. the third body. Ding [START_REF] Ding | A finite element based approach to simulating the effects of debris on fretting wear[END_REF] tried to take this into account in a new finite element model (without representing this third body), but the wear tracks that can be simulated are shallow. Semi-analytical approaches are also used [START_REF] Gallego | A fast and efficient contact algorithm for fretting problems applied to fretting modes I, II and III[END_REF][START_REF] Gallego | Multiscale computation of fretting wear at the blade/disk interface[END_REF], which allow reducing the computation time and simulate the full system. The reader can refer to the article [START_REF] Renouf | Numerical tribology of a dry contact[END_REF] for more information on the numerous methods used in finite elements.

Although finite element models allow modeling the stress and displacement fields of bodies in contact, they do not allow considering the discontinuities of the material. The discrete element approaches used in tribology have their origins in the numerical modeling of discontinua of Cundall and Hart [START_REF] Cundall | Numerical modelling of dicontinua[END_REF]. These simulations, while quite varied, are based on three key principles:

• The equations of motion of mechanics;

• The detection of contacts between particles;

• The calculation of the forces between the particles. Discrete element models (DEM) are used in tribology to model different aspects of contact. For example, DEM can be used to model the contact and the degradation of materials between the wheel of a train and the rail [START_REF] Chapteuil | Influence of copper/graphite properties on the tribological and electrical behavior of copper-graphite third body layer[END_REF]. DEM simulations are also used to simulate the flow of the third body-as well as its degradation-in different contexts: thermo-mechanical response [START_REF] Rivière | Thermo-mechanical investigations of a tribological interface[END_REF], two-phase materials [START_REF] Champagne | Modeling wear for heterogeneous bi-phasic materials using discrete elements approach[END_REF], or electrical contacts [START_REF] Renouf | Coupling electrical and mechanical effects in discrete element simulations[END_REF][START_REF] Descartes | A new mechanical-electrical approach to the wheel-rail contact[END_REF].

Some authors tried to combine discrete elements with finite elements to combine the advantages of both methods [START_REF] Renouf | First-body versus third-nody: dialogue between an experiment and a combined discrete and finite element approach[END_REF][START_REF] Li | Fretting damage modeling of liner-bearing interaction by combined finite elementdiscrete element method[END_REF][START_REF] Leonard | Third body modeling in fretting using the combined finitediscrete element method[END_REF]. In particular, the particles of a discrete model are a priori rigid, which is not the case when combining these two methods. The field of deformable granular materials is of growing interest, both experimentally [START_REF] Dijksman | Refractive index matched scanning and detection of soft particles[END_REF][START_REF] Vu | Soft-grain compression: beyond the jamming point[END_REF] and numerically [START_REF] Cantor | Compaction model for highly deformable particle assemblies[END_REF][START_REF] Nguyen | Compaction of granular materials composed of deformable particles[END_REF][START_REF] Harthong | Contact impingement in packings of elastic-plastic spheres, application to powder compaction[END_REF][START_REF] Favier De Coulomb | Rheology of granular flows across the transition from soft to rigid particles[END_REF], and extends well beyond the purely tribological applications.

Finally, the last type of approach concerns nanometric models. Several techniques exist, based on the displacement of atoms. The Monte Carlo method and the molecular dynamics method are the most widespread. These simulations are based on atomic interaction laws derived from quantum mechanics [START_REF] Doucet | Computer-Aided Molecular Design: Theory and Applications[END_REF]. Some authors focus on modeling adhesive wear and the formation of agglomerates in contact [START_REF] Aghababaei | Critical length scale controls adhesive wear mechanisms[END_REF][START_REF] Molinari | Adhesive wear mechanisms uncovered by atomistic simulations[END_REF]. These simulations are informative but timeconsuming, and only allow representing a tiny fraction of the contact. For more information on these simulations, the reader may refer to [START_REF] Robbins | Computer simulations of friction, lubrication, and wear[END_REF].

Methods

Numerical model

This work is largely based on a previous work [START_REF] Mollon | Solid flow regimes within dry sliding contacts[END_REF], where an extensive simulation campaign was performed, with various mechanical parameters. In what follows, we revisit these stored numerical results with dedicated techniques to extract relevant new kinematic data. It is therefore worth providing some general information about the model. The reader interested in more details may refer to [START_REF] Mollon | Solid flow regimes within dry sliding contacts[END_REF]. In that paper, the author describes a simulation of a flow of deformable particles in a dry sliding interface. These simulations are performed on the open-source code MELODY, a software developed in LaMCoS, INSA Lyon, France [START_REF] Mollon | A multibody meshfree strategy for the simulation of highly deformable granular materials[END_REF][START_REF] Mollon | A unified numerical framework for rigid and compliant granular materials[END_REF]. This software generalizes the conventional Discrete Element Method (DEM) and allows simulating complex granular materials, with either rigid or highly deformable grains with arbitrary shapes, within the framework of large deformation hyperelasticity in 2D. The versatility of this software allows it to be used in several fields, for example in mechanical engineering [START_REF] Zhang | Significance of third body rheology in friction at a dry sliding interface observed by a multibody meshfree model: influence of cohesion between particles[END_REF][START_REF] Quacquarelli | A dual numerical-experimental approach for modeling wear of Diamond Impregnated Tools[END_REF] or geomechanics [START_REF] Casas | DEM analyses of cemented granular fault gouges at the onset of seismic sliding: peak strength, development of shear zones and kinematics[END_REF].

Principles

Figure 2 shows a sketch of the numerical model, which is purely dimensionless. Two rough rigid walls (length L = 100 ) are placed on either side of a collection of 2000 deformable particles with a unit average diameter, to represent two bodies rubbing against each other with a third body. The two first bodies have a self-affine roughness, with a R a equal to 1 space unit. They are rigid and non-degradable. Their cohesion with particles is very high, in order to prevent wall-slip effects, and to focus on velocity accommodation within the third body only.

The deformable particles follow a hyperelastic law. This makes it possible to apply elasticity to bodies that undergo large deformations and rotations. The neo-Hookean law considers two main parameters, which are Young's modulus E and Poisson's ratio . The Poisson's ratio is fixed at 0.495 to impose a quasi-incompressibility. Finally, the density is equal to 1 (dimensionless model). To make the simulation stable, damping is added via the Rayleigh parameter: the damping matrix is proportional to the stiffness matrix via this parameter. The result is a viscoelastic-like behavior of the Kelvin-Voigt type. The contact law governing the interactions between the particles and with the walls is governed by a cohesive frictionless model with only one parameter c , which is the strength per length unit of contact between two surfaces. By using deformable particles that can move freely in space, the model combines the laws of dynamics with those of continuum mechanics.

A gap (noted D ) between the two first bodies is imposed. Finally, a relative tangential velocity V = 1 is imposed on the two bodies: the upper body goes rightwards at V∕2 , and the lower one, leftwards at -V∕2 . This generates a shear within the collection of particles that will accommodate velocities while transmitting loads throughout the simulation, like the third body within a tribological contact. Periodic boundary conditions are applied to the extremities of the simulation. For each simulation, the simulated time corresponds to three sheared spatial periods, i.e. a time of 300-time units (called Δt and equal to 1 in this dimensionless context). To ensure that the flow regime of the third body is well established, only the last 200Δt will be considered for quantitative interpretation.

The normal and tangential forces experienced by the upper and lower bodies are measured over time and noted F N (t) and F T (t) respectively. Their time average is noted F N et F T . The stresses, which are the ratio of the forces to the length L (2D simulation), are noted N and T . The ratio of these two quantities provides an average coefficient of friction at the scale of the whole interface, as an output quantity of the simulation. This approach is a good conceptualization of the experimental reality where the friction coefficient of a given interface is a measured quantity that emerges from a large number of local phenomena and events. The elastic behavior of each grain can certainly be seen as a limitation when applied to plastic and ductile materials (e.g. metals), but constitutes an improvement with respect to perfectly rigid grains. Besides, irreversible deformations occur within the flow because of the cohesive character of the contact laws. It means that plasticity is present in the simulations, albeit at the scale of the flow instead of that of the grains.

Parametric space

Three dimensionless parameters are studied and lead to various flows of the third body.

These three parameters are a normalized stiffness Ẽ (between 0.15 and 1.2), a normalized cohe- sion c (between -1 and 1.4) and normalized viscosity ̃𝛼 (between -3.3 and -1.5). The first two parameters are normalized by the average normal stress experienced by the first bodies, n . The third parameter is normalized by the mean strain rate in the interface, V∕D . The normalized stiff- ness characterizes the resistance to deformation of particles under stress. The normalized cohesion quantifies the cohesive force between the particles in the system. This cohesion therefore varies roughly between one order of magnitude below the normal stress (negligible cohesion) and one order of magnitude above it (dominating cohesion). Finally, the normalized viscosity characterizes the response of the grains to their rate of deformation. As demonstrated in [START_REF] Mollon | Solid flow regimes within dry sliding contacts[END_REF], the chosen range for ̃𝛼 is consistent with the solid viscosity of steel (for example) under typical tribological conditions.

Consequences of stiffness, cohesion, and viscosity of third body

The coefficient of friction, noted , represents the ratio of the measured tangential force to the normal force. This is an output value, not an input value. It characterizes the shear strength of this confined soft flowing granular material.

The evolution of the friction coefficient as a function of the values of Ẽ , c and ̃𝛼 can be seen in Fig. 3. The coef- ficient of friction has a wide range, from 0.1 to over 1.0, but remains in the typical range of observed values for dry sliding. The increase in the ̃𝛼 parameter leads to a general increase in the coefficient of friction. The two maps of Fig. 3 show the coefficient of friction at imposed ̃𝛼 . In particu- lar, when ̃𝛼 ≈-1.8 , a maximum coefficient of friction is observed, at Ẽ = 0.7 and c = 0.7 (that is, for a Young modu- lus of the grains and an interparticle cohesion roughly one order of magnitude larger than the applied normal stress).

These three parameters greatly influence the flow and the microstructures of the third body. Figure 4 shows different kinematic patterns of the third body as a function of the parameters Ẽ and c , at imposed ̃𝛼 =-1.8 . For example, when Ẽ and c are low, the observed behavior is close to

(4) = F T F N picture A.
The particles deform a lot and fill the pore space. They adopt a structured pattern with a preferred horizontal direction which promotes a perfectly laminar flow: because of the low cohesion, each layer of particles slides over the others. On the other hand, in the case of B, the particles deform little. The observed regime can be described as granular. Cases D and F are also particularly interesting. In these cases, the particles have sufficient cohesion to form agglomerates, similar to those seen in Fig. 1. These third body agglomerates flow, shear and roll into the interface. The other cases (C, E, H, and I) are intermediate to those described above. These simulations are quantitatively consistent with the experimental variety of third body visual aspects [49].

Coherence

In some simulations, agglomerates are observed. Grains form various sizes of cohesive clusters, in which the particles move in a motion similar to that of a rigid body. An example can be seen in Fig. 5. The agglomerates are difficult to observe in a succession of still images but are clearly visible in animations.

The characterization of the size of these agglomerates would make it possible to analyze the consequences of the rheology of the third body particles on the first bodies. Tools developed to characterize this size of agglomerates, and derived from fluid mechanics are documented and adapted hereafter.

The agglomerate identification function is derived from vortex identification functions. These mathematical descriptors were proposed in the context of vortex studies in fluid mechanics and are called Γ 1 and Γ 2 [50, 51]. These functions compare the observed motion of points near a central point with that of particles forming an ideal vortex. The function proposed here, called "coherence", is based on the same concept but compared with the motion of an ideal rigid body. To do this at a given time step, a velocity field is first interpolated on a regular Eulerian grid from the motion of the particles, as can be seen in Fig. 6. The interpolation used is cubic, which allows C 2 continuity. The velocity is interpolated with respect to the particle center-of-mass velocity, because the deformation-related velocity field in particles is small compared to the rigid body velocity of these particles, as can be seen in Fig. 2. The time step used to extract grains positions and to compute velocities is equal to one time unit, and corresponds to an average strain level of 0.05 between two extractions. This is a good compromise between too small (producing a large amount of redundant data) and too large (overlooking complex trajectories of the grains) time steps. The velocity field is discretized into a square mesh, with a spatial step of Δx =Δy = 1 , and each point of this mesh has a velocity ⃗ V . We then define a local quantity for each point, called "coherence", which quantifies to what extent the velocity field in the neighborhood of this point corresponds to a rigid-body motion. Coherence at a given point P is therefore calculated by comparing, at each point P , the velocity of each point M located in a disk of radius R and center P , with the velocity that an ideal rigid body of center P and radius R would have. The complete formula is as follows:

where PM ≤ R , V circ is an ideal circumferential relative velocity of point M with respect to point P in a rigid-body rotation (the angular velocity is arbitrarily chosen to be equal to 1), and V ave is the average velocity of all points M within the radius R (cf. Fig. 7). By construction, -1 ≤ C < 1 . A value of C = 1 means the particles rotate in a coordinated rigid-body motion in the positive direction, and vice versa for C =-1 . The radius R of the interrogation circle is a user- defined parameter. It seems that a radius of two to three times the typical size of a particle allows the detection of structures properly while keeping the computation time reasonable.

(5)

C(x, y, R)= ∫ M∈S V -V ave ⋅ V circ dS √ ∫ M∈S V -V ave 2 dS ⋅ √ ∫ M∈S V circ 2 dS
By repeating the operation on the whole simulation, it is thus possible to determine a coherence field, as can be seen in Fig. 8A, computed with the velocity field extracted from Fig. 5. Red zones correspond to a coherence close to 1 , and blue ones to -1 . Expectedly, the coherence map is essentially negative because the natural direction of rotation due to the shearing of the interface is negative in the present sign convention. By applying a filter to select only the areas where |C| > 0.8 , agglomerates are highlighted (Fig. 8B). It is finally possible to stack these detected areas to observe their evolution over time (Fig. 8C). Independent structures are easily identified and labeled based on connectivity relationships in the space-time domain.

Each detected agglomerate is characterized by its coherence. An additional relevant property to be calculated is the average angular velocity, noted . Each agglomerate detected at an instant T can be associated with a certain number of grains belonging to it. The average angular velocity is calculated by computing the average of the circumferential component of the velocity of each grain to the barycenter of the agglomerate.

With R i the distance between each point in the structure to the barycenter and n the number of points in. To consider the difference in distance between the first two bodies according to the simulations, this quantity is normalized by the largescale shearing of the third body layer, V∕D , and noted ̃𝜔 .

In this case, the V∕D shear rate is negative. Structures rotating in a positive direction will therefore have a positive value of ̃ . Normalization makes it easy to compare the angular velocity of the structures with the shear rate. If ̃𝜔 = 1 , the structure rotates as a disk in a rolling movement without sliding in the contact. If ̃𝜔 > 1 , it rotates faster, and if 0 < ̃𝜔 < 1 , the structure rotates slower. Finally, if ̃𝜔 < 0 , the structure rotates in the opposite direction to the shear.

Results

75 simulations were performed with different values of Ẽ , c and ̃𝛼 in [START_REF] Mollon | Solid flow regimes within dry sliding contacts[END_REF], and we revisit their stored numerical results. Coherence is applied to this series of simulations used to evaluate the coefficient of friction as a function of the three parameters identified above. Seven notable simulations are retained and are the same as those used in the previous study. They are noted from A to G and located in the parametric space in Fig. 4.

(6) (x, y, R)= ∑ V circ R i � n (7) ̃𝜔 =
𝜔 V∕D Fig. 7 Coherence calculation scheme

Stacked coherence

In Fig. 9, the coherence is plotted over time and stacked, following the methodology shown in Fig. 8C. In case A, the particles are soft and hardly cohesive. There are very few agglomerates detected, and they all rotate in the negative direction. Because of the shear generated by the two upper and lower bodies, this is naturally the preferred direction of rotation, while positive rotations correspond to agglomerates which rotates against the main shear direction. The lifespan of the structures is very short.

Cases B, C, and D follow a path of decreasing stiffness. C is located at the maximum of the friction coefficient. The other parameters do not vary very much. First of all, for case B, structures are not very spatially extensive but persist for a fairly long time. No observed structure rotates in the positive direction. For case C, on the other hand, the structures, although numerous, are not very extensive spatially and temporally. Some structures rotating in the positive direction are to be noted, unlike case A. Finally, in case D, the negatively rotating structures are massive, both spatially and temporally. Structures that rotate in the positive direction are also observed, and their spatial extent is notable. We thus observe for these three cases, with decreasing Ẽ , first many struc- tures, then a drastic reduction in their size, and finally the appearance of large, numerous, and persistent agglomerates.

Simulations E, C, and F follow an increasing cohesion path. Case E is similar to case A, except that there are more structures, and they have a slightly larger spatial extent. Case F is very interesting and shows the impact of an important cohesion. The structures are very large, whatever their direction of rotation. In these three cases, the increase in cohesion results in an increase in structures size. Cases E and C are similar, while regime F has very large structures, but is seemingly shorter-lived than in case D.

Case G, visible in Fig. 10, is similar to case C, except that the viscosity is much lower. The alpha viscosity has an important role in whether or not structures are formed since the detected structures have a longer lifetime and a larger spatial extent.

Size of structures

Figure 11 shows the proportion of grains belonging to a structure compared to all the grains in the simulation. This of grains is very dependent on the cohesion c , and to a lesser extent, on the viscosity ̃𝛼 . On the other hand, the stiffness Ẽ seems to have very little influence. The maximum grain content in a structure, 11%, is rather low. This is due to the choice of a rather high detection threshold ( C > 0.8 ). The data should be analyzed relatively rather than absolutely. The proportion of grains in the structures is an important parameter since it gives information on the extent of the agglomeration phenomenon.

Angular velocity of structures

The probability density function (PDF) is used to extract statistical data from all the simulations, including the angular velocity of the agglomerates. These PDF are in terms of equivalent grains, i.e. they represent the density of probability for a given grain to belong to a structure with a certain diameter or a certain angular velocity. They are visible in Fig. 12, left column. The seven cases analyzed above have very different profiles.

Path B (purple) → C (medium blue) → D (green), which follows a decreasing stiffness, is interesting because we observe that case C (maximum peak of the coefficient of friction) corresponds to a minimum angular velocity and a minimum of grains in the structures. For path E (dark blue) → C (medium blue) → F (light blue), which follows an increasing cohesion, the case C is intermediate between E and F. Case G does not have very different characteristics from case C, which has the same stiffness and cohesion, but not the same viscosity. Case A, in orange, which has very few structures, shows only structures rotating in the shear direction, with very high dispersion. Thus, many structures rotate nearly twice as fast as the shear rate.

To each of the 75 simulations can be attributed an angular velocity mode ̃𝜔 corresponding to the peak of the angular velocity of the PDF function. Some correlations are interesting to observe, for example, between Ẽ and ̃𝜔 in Fig. 13A. No correlation (Pearson's coefficient of 0.22, Pearson's coefficient is relevant to quantify the quality of linear regression) is noted between ̃𝜔 and Ẽ for the whole simulations, but a correlation is visible when the data are binned according to the value of the normalized viscosity. The linear relationship is much stronger at a low viscosity (Pearson's coefficient of 0.51 for ̃𝛼 =-2.2 and 0.55 for ̃𝛼 =-3.2 ); at ̃𝛼 =-1.8 , the correlation is weaker. Another remarkable correlation, visible in Fig. 13B, is that between the average angular velocity of the structures ̃𝜔 and the solid fraction S F . The solid fraction is the ratio of the surface occupied by the particles in the simulation to the surface of the space that contains them. A solid fraction close to 1 indicates that the grains are very tightly packed and fill all the pore space, while a value close to 0 indicates the opposite.

Diameter of structures

One way to quantify the size of the structures is to use a cylindrical representation. For this purpose, each structure is associated with a lifetime L and a diameter D. This diameter D is based on the average area of the structure in each time step along the lifetime. It is thus possible to approximately represent each structure as a cylinder whose height is the lifetime and whose diameter characterizes the average spatial extent over time. An example of cylindrical representation can be found in Fig. 14.

In the same way as for the normalized angular velocity, statistical data of the structure diameters can be extracted from the studied simulations. The probability density function and the cumulative distribution function of the diameter of structures can be seen in Fig. 12, right column. In path B → C → D, case C corresponds to a minimum size of agglomerates. In the same way, as for the angular velocities, case C corresponds to an intermediate step in the E → C → F paths.

An interesting relationship can be observed between the cohesion and the characteristic (i.e. the mode of the PDF) Fig. 13 A ̃ as a function of normalized rigidity Ẽ ; B ̃ as a function of solid fraction S F . The size of the dots is relative to the proportion of grains in agglomerates in the simulation (in %) Fig. 14 Cylindrical representation of detected structures for case F (cf. Fig. 9). The structures rotating in the positive direction are in red, and in blue for those rotating in the negative direction. The intensity of the color depends on the angular velocity diameter of the structures (cf. Fig. 15). The Spearman's correlation coefficient (relevant for non-linear correlations) is quite high, over 0.67. This indicates that the larger the cohesion, the larger the diameter of the structures, which agrees with the observations that can be made in Fig. 9. The relation between these quantities is non-linear, with a sharp increase of the characteristic diameter of the agglomerates when the normalized cohesion goes beyond 0.5.

Coefficient of friction

A correlation can be found between the normalized angular velocity structures ̃𝜔 , and the coefficient of friction and can be visualized in Fig. 16. The Pearson correlation coefficient between these two data is -0.70. When a large friction coefficient is observed, the average angular velocity of the field is very low. In contrast, the friction coefficient is found to be very low when the normalized angular velocity gets closer to 1, i.e. when the agglomerates rotate following the natural shear rate of the flow.

As described in [START_REF] Mollon | Solid flow regimes within dry sliding contacts[END_REF], based on energetic arguments, coefficient of friction can be expressed as a sum of two partial friction coefficients: a "surface-related" coefficient of friction S , and a "bulk-related" coefficient of friction B . The former describes the part of the energy dissipated during the creation and destruction of surfaces in the granular sample during shear, and the latter the energy dissipated by inelastic deformations within the deformable particles. These components are particularly complex and difficult to correlate with other parameters. However, two significant correlations can be observed. In Fig. 17A, a link can be made between the normalized angular velocity of structures ̃𝜔 (cf. Sect. 3.2) with the bulk-related coefficient of friction, with a Pearson coefficient of -0.58. As the structures rotate more, there is less energy loss through dissipation in the bodies.

In Fig. 17B, a clear link can be made between the size of the structures (cf. Sect. 3.4) and the bulk-related coefficient of friction (Spearman coefficient of -0.83). This relationship is monotonic but non-linear. The higher the coefficient of friction associated with the interparticle energy dissipations, the smaller structures there are. In contrast, the appearance of large agglomerates (i.e. with D > 4 ) cor- relates well with a large drop of the surface-related energy dissipation in the system.

Discussion and conclusion

Based on the quantitative results presented in the previous section, we can draw a general picture of the causes and consequences of the agglomeration phenomena in soft cohesive granular flows. It appears that the existence of persistent rigid-like agglomerates in this class of flows is promoted by the interparticle cohesion and to a lesser extent by the viscous response to shear deformation (cf. Fig. 11). Agglomeration is indeed negligible when the quantities c and ̃𝛼 are low. The size of the agglomerates is mostly related to cohesion (cf. Fig. 15), which is simply explained by the fact that more cohesive particles tend to pull more surrounding particles in the same motion. The normalized cohesion value c = 0.5 seems to be a threshold above which very large structures can develop. Viscosity also plays a role in the creation of large structures (cf. Fig. 11).

The stiffness and the viscosity parameters play an important role in the angular velocity of the structures (cf. Fig. 13A). The softer and less viscous the grains are, the more they tend to form structures with an angular velocity corresponding to the natural shear rate of the interface (although a few outlying cases are noticed). When the stiffness is important, the structures are more reluctant to follow the natural shear-related rotation. The solid fraction also plays a role in the angular velocity of the structures (cf. Fig. 13B). When the solid fraction is high, the grains are closer to each other, their contacting areas are larger, and their motions are more coordinated, which seems to increase the angular velocity of agglomerates. Since a high stiffness and a low viscosity (to a lesser extent) decrease the solid fraction, we can speculate that this decrease in solid fraction provides less incentive for structures to rotate. Hence, a simplistic view of this class of granular flows could be that cohesion and viscosity favor the formation and the growth of agglomerates, and that stiffness and viscosity restrict their natural rotation within the sheared flow. The proposed results hence provide a contribution to the growing topic of the rheology of cohesive granular flows [START_REF] Macaulay | Viscosity of cohesive granular flows[END_REF][START_REF] Mandal | Rheology of cohesive granular media: shear banding, hysteresis, and nonlocal effects[END_REF], with the added complexity of a cohesive strength increasing with the contact area between the grains (and hence with their deformability and the applied stress).

A link can be made with the surface-related coefficient of friction (cf. Fig. 17B). At a diameter greater than 4, the value of the surface-related coefficient of friction drops and becomes very low. This could be due to a geometrical effect linked to a specific surface limit beyond which the behavior changes radically. This value is approximate and could be dependent on the parameters chosen in the simulations. It is at the same value that a behavior change is observed between the cohesion and the average diameter of structures. This decrease of the surface-related dissipation can be explained by the fact that the agglomerates become larger. The specific surface area thus decreases in the sheared sample and the contact area mobilized between the agglomerates diminishes. The energy dissipation therefore no longer occurs at the interfaces between the agglomerates, but within them, via the inelastic component of their deformation. Likewise, as structures rotate faster (i.e. closer and closer to the natural angular velocity associated with rolling in the sheared flow) they tend to endure less deformation, resulting in a decrease in the bulk-related coefficient of friction (cf. Fig. 17A).

Finally, the friction coefficient strongly depends on the angular velocity of the agglomerates (cf. Fig. 16). It can Fig. 17 A Bulk-related coefficient of friction as a function of mean normalized angular velocity; B Surface-related coefficient of friction in function of mean diameter be assumed that the agglomerates roll in the contact, and decrease the coefficient of friction, like balls in bearings. This link between high cohesion, the generation of rolling agglomerates, and a limitation of the resulting friction coefficient at the scale of the interface is consistent with the hypothesis reported in a previous study [START_REF] Zhang | Significance of third body rheology in friction at a dry sliding interface observed by a multibody meshfree model: influence of cohesion between particles[END_REF], but is now supported by more quantitative data. The accommodation of the relative velocity of the two walls is more likely to occur via rolling at low stiffness and high cohesion, while it is more likely to occur via particle deformation when the cohesion is low.

The methodological framework proposed in this paper for the detection and quantification of agglomerated structures is still imperfect. A blind spot exists for structures with zero angular velocities, although it seems minimal. The presence of large pore spaces is also ignored by the interpolation algorithm, and this could be corrected in future work. The main limitation, however, seems to be the rather short duration of the simulations, which are imposed by the large computational cost of the multibody meshfree approach but limit the statistical accuracy of the detection and characterization of the agglomerates. Much longer simulations could be made possible in future studies by the use of the Soft Discrete Element Method [START_REF] Mollon | The soft discrete element method[END_REF]. The thickness of the sheared layer, which provides an upper bound for the agglomerates diameters, might also play an important role and should be investigated in future studies.

From a more applicative perspective, the tool developed here will be applied to fretting simulations in a tribological context. By applying this approach to realistic mechanical contacts, we expect to establish a link between the third body rheological properties, the formation and characteristics of agglomerates, the stress fluctuation patterns undergone by the surfaces bounding the third body flow, and finally their damage and wear.

Fig. 1

 1 Fig. 1 SEM views of different areas of mechanical contact after friction tests, with various third body accommodation regimes

( 1 )Fig. 2

 12 Fig. 2 Sketch of the numerical model with two rigid first-bodies and a third-body represented by deformable grains. Colors are arbitrary and only serve to distinguish grains on the top and the left figure. Velocity field magnitude is plotted in the lower-right corner

Fig. 3

 3 Fig. 3 Coefficient of friction in function of Ẽ and c , for ̃𝛼 =-1.8 and ̃𝛼 =-3.2

Fig. 4 Fig. 5 Fig. 6

 456 Fig. 4 Microstructure of third body as a function of Ẽ and c (from [43])

Fig. 8 A

 8 Fig. 8 A Coherence field extracted from the velocity field (Fig. B Filtered coherence field at the same time-step at |C| > 0.8 ; C Stacked filtered coherence over time

Fig. 9 A

 9 Fig. 9 A-F Stacked filtered coherence over time for six characteristic simulations. The structures in red rotate in the positive direction, and vice versa for the blue ones

Fig. 10

 10 Fig. 10 Stacked filtered coherence over time for G case. The structures in red rotate in the positive direction, and vice versa for the blue ones

Fig. 12

 12 Fig. 12 Probability density function of structures for various cases of normalized angular velocity of structures ̃𝜔 and diameter of structures D

Fig. 15 Fig. 16

 1516 Fig. 15 Average diameter of the structures as a function of the cohesion parameter

Safran Aircraft Engines, Moissy-Cramayel, France

Acknowledgements

The authors gratefully acknowledge Safran Aircraft Engines and the French National Research and Technology Agency (ANRT) for financially supporting this research project. The authors acknowledge that this study contains original material. Its publication has been approved tacitly by the responsible authorities at the institutes where the work has been carried out.

Declarations

Conflict of interest The authors declare that they have no conflict of interest.