Force fluctuations in a pushed granular material

We analyse the time evolution of stresses in a noncohesive granular bed cut by the motion of a tool. Our numerical simulations both in two and three dimensions reveal large fluctuations of the force acting on the tool. These fluctuations have a decreasing exponential distribution. We find that, in spite of fluctuations, the mean force is well fitted by an analytical form obtained from a limit analysis with Coulomb’s failure criterion.


Introduction
In several technological domains involving granular materials, such as soil mechanics and silo engineering, a popular approach to the numerical prediction of the overall response of a granular mass to external actions relies on modelling it as an elasto-plastic continuum with a failure criterion. Finite element calculations are then applied to the corresponding boundary value problems [1][2][3]. In recent years, however, the numerical simulations of the dynamics of large collections of solid particles have become efficient enough to allow for the solution of such problems of overall response on the basis of grain-to-grain interactions. These methods, known as discrete element methods [4][5][6], yield at the same time the macroscopic description of the behaviour of the investigated granular mass and some microscopic details about local phenomena. An obvious advantage of this approach is that it needs no other physical input than the interaction laws between the particles and, in this respect, it provides an excellent way of testing the validity of macroscopic models and their underlying assumptions. This study is part of a research work on mine removal, carried out within a co-operative program devoted to surface soils (PEA Sols de Surface) developed in France under the heading of ETAS (Etablissement Technique d'Angers, DGA). We acknowledge its financial support (contract N. 9656906).
This paper is concerned with the numerical investigation of the action of a simple rigid "tool" upon a granular material. In 2D numerical experiments, the tool is partially immersed in a granular bed composed of irregular polygonal particles initially in equilibrium under gravity inside a rectangular box. Then the tool is moved horizontally with constant speed. In 3D numerical experiments, the tool consists of a rigid assembly of small overlapping spherical particles, reproducing the shape of a plane rectangular blade.
We are interested mainly in the total force acting on the tool. Several authors have reported on large temporal stress fluctuations in silos [7,8] or in a sheared granular material contained in a box [9,10]. The physical origin of these large fluctuations still remains obscure. They might simply be a reminiscence of the highly heterogeneous distribution of contact forces in granular media [11][12][13] or result from kinematic constraints due to particular boundary conditions. In our 2D and 3D simulations, we similarly found large force fluctuations while the external forcing seems far less stringent than in the situations reported above. In fact, our tools interact with the granular mass through a layer whose thickness does not exceed a few particle diameters. We will show that a stretched exponential form provides an excellent fit for the statistical distribution of these fluctuations. Another significant result of this study is that, in spite of the fluctuations, the mean force is very well fitted by the analytical expression obtained by standard limit analysis calculations based on the Coulomb failure criterion.
The simulations reported here were part of a technological research project devoted to the mechanical behaviour of surface soils, with a view to the design of soil-cutting tools such as those affixed to vehicles intended to remove mines [14]. In this paper, we will focus on the evolution of forces. The main results will be presented after a short description of the simulation method and of the physical parameters of our simulated systems. Mean forces and force fluctuations will be discussed in two sections. Several aspects of the problem such as the influence of system size and cutting speed will be addressed as well.

Simulation method
For these simulations, we used the Contact Dynamics (CD) method. For a detailed description of this numerical technique the reader may refer to several papers published in recent years [15][16][17][18]. The method consists of an implicit time-stepping resolution of the dynamical equations of a collection of rigid particles interacting at their contact points. The condition of impenetrability and the dry friction law of Coulomb are implemented without resorting to the smoothing schemes used in other discrete element algorithms [5]. Contact Dynamics has been successfully applied to study both dynamical phenomena, such as size segregation [19], and quasi-static deformations of granular media [20,21].
In 2D simulations, we used circular or polygonal particles and, in 3D, spherical ones. Since in 2D the main trends of the results were found nearly independent of particle shape, the 2D simulations presented in this paper all rely on particles of randomly generated convex polygonal shapes with at most 8 edges. Once a polygonal shape is generated, the program determines its mass-centre (corresponding to uniform mass-density) and the control disk i.e. the smallest disk centered at this point containing the polygon. In the experiments presented here, the control radii exhibit rather small dispersion (see Table 1). All samples are prepared by randomly sprinkling grains in the presence of gravity into boxes. Before starting a numerical experiment, one leaves the deposited mass to relax until the competition of gravity and intergranular Coulomb friction leads to a state of sufficiently small kinetic energy. In 2D, we have four samples of 1110, 2039, 4598 and 10437 grains (which occupy the same volume) to which we will refer as samples A, B, C, and D, respectively. In 3D, in order to save computation time, we have assumed that the flow, as well as the tool during its whole motion, admit a plane of symmetry. The latter is realized as a frictionless rigid boundary. Moreover in the 3D sample, to which we refer as sample E, a layer of grains is glued to the internal walls of the container (see Fig. 1). This artificial roughness is meant to prevent horizontal recirculation of the particles due to the finite size of the box. The control diameters and void ratio of all the numerical samples are given in Table 1. The interparticle coefficient of friction equals 0.3. Walls, as well as the glued grains, are treated as frictionless.
Tools are rigid assemblies of polygonal or spherical particles. A tool may have different shapes and different inclinations to the vertical. The simplest tool we used in 2D was a vertical beam composed of a rectangle completed with an octagonal tip in order to make the cutting edge smoother. Their coefficient of friction with the free particles composing the samples was 0.3. In each numerical experiment, the tool was first driven vertically into the sample with a constant speed until the desired depth was reached. Then, it was moved with a constant horizontal speed v along the granular bed. Most simulations were performed with v = 2 cm/s. A few of them, however, were performed with cutting speeds of 50 and 100 cm/s. Figure 2 shows three stages of a typical experiment with v = 2 cm/s. The time-step was 0.1 ms in all 2D simulations. In 3D, since we are mainly interested in the pushing phase, we suppressed the initial penetration phase. The blade was just set on the outside of the box, as near as possible to the granular sample, the program not taking into account any visibility between tools and pavement grains. The tool was then dragged horizontally with a constant speed v = 2 cm/s with a time step equal to 0.2 ms. We have chosen the same coefficient of friction between tool and free grains as in 2D, i.e. 0.3. Figure 3 shows one stage of a such 3D computation, in which we can observe the lateral flow of grains.

3
The mean force Figure 4(a) shows the typical evolution of the horizontal and vertical components of the force F acting on the tool as a function of time in a 2D simulation. The simulation is stopped when the tool comes too close to the right wall. The whole episode of horizontal displacement of the tool takes 100000 time-steps. Two features are evident in the force signal. The force increases at long time scales, due to the pile growing in front of the tool, but it is affected by large fluctuations at shorter time scales. Figure 4(b) shows the force variation in a 3D simulation. The fluctuations here seem less important than in the 2D case. In 3D, the pile in front of the tool reaches a quasi-stationary size at about t = 1 s, due to the possibility of lateral flow and the force remains almost constant after that. It further increases as the tool comes too close to the right wall. In Fig. 5 we have plotted the horizontal component of the force exerted by the tool on the grains, averaged over intervals of 2000 time steps, corresponding to 0.2 s, for the samples A to D. These four 2D samples which differ in coarseness have the same overall geometry and the penetration depth of the tool is the same. The mean force is almost the same for all four samples at the respective stages of the evolution of the system. We observe also low frequency fluctuations of the average force. What does the mean force represent for this system? As we shall see below, large deviations from the mean force are exponen- tially rare. But, we checked that the mean force is nearly independent of an upper limit over the forces taken into account for its calculation. Thus, the mean force is an integrated quantity of well-represented forces over a given time interval.
In order to fit some analytical expressions to the evolution of the mean horizontal force F , in what follows we consider two different models. The first one is Coulomb's method of wedges [22] which yields the following expression for the force exerted by a pushed granular bed of horizontal top surface against a wall: where γ is the bulk density, h is the depth below the top surface, φ is the internal angle of friction, and k = µ w / tan φ with µ w denoting the global coefficient of friction between the tool and the granular material. This expres- where the granular medium resists the imposed displacement. In the opposite case, this is the tool that resists the action of the medium. In this active state, the denominator in the above expression should be replaced by (1 + √ 1 + k sin φ) 2 . Of course, this expression cannot directly be applied to our problem since the top surface in front of the tool is inclined to the horizontal. Although the shape of the pile is not exactly that of a pile at its angle of repose, it is useful, as a first approximation, to simplify the geometry of the system to the one displayed in Fig. 6. We then apply the expression (1) to an "effective depth" h given by h = 2 3 h 1 + 1 3 h 2 . This corresponds to a top line passing through the center of gravity of the pile. In other words, the upper half of the volume of the pile is redistributed on the top surface. This first choice is not really justified, but seems reasonable in view of the following results. On Fig. 7 we have shown the mean force calculated in simulations and the fit obtained by this simple model as a function of h 2 both for active and passive states. We see that the mean simulated force stays very close to the passive force. For this fit φ is evaluated from the angle of repose of the half-pile left behind the tool as it moves forward; one finds φ = 27 • . As for µ w , it is estimated from the direction θ of the resultant force on the tool. This angle fluctuates in time, but its most probable  value is θ = 17 • with horizontal, which corresponds to µ w = 0.3.
Let us now consider a rigid-plastic calculation of the passive force for the geometry represented in Fig. 8. This has been done by Perumpral et al. [23] for 3D systems and for a tool inclined at an angle β to the horizontal and with an angle of friction δ and for the Coulomb failure criterion with cohesion. This model assumes that the deformation is localized on a plane passing through the bottom of the tool and making an angle ρ with the horizontal. This model yields the following simplified expression for the passive limit of horizontal force in a 2D cohesionless material: where W is the total weight of the mobilized domain, i.e. above the sliding surface, given by Fig. 8. Geometry of the mobilized domain considered by Perumpral et al. according to [23] where h 1 and h 2 have the same meanings as in Fig. 6. In the 3D case, the expression of the drag force is given by the following expression: in which the friction force SF 2 due to one of the lateral part of the pad is taken into account: In this case, the total weight W is given by where B is the total width of the blade. Our simulations show that most of the deformation occurs, as expected, along a line passing through the bottom of the tool as it can be seen in Fig. 9. Each particle has been coloured with a level of grey corresponding to the local shear rate of the flow in the vicinity of this particle. This shear rate is estimated by fitting a velocity gradient tensor to the distribution of the velocities of the centers of the neighbouring particles. One could use the location of this major shear zone in order to estimate ρ. Another possible procedure is to minimize F in Eq. (2) with respect to ρ as suggested by Perumpral et al. [23]. This gives ρ = 20.7 • independently of h 2 in the 2D case. This is a plausible estimate of ρ in view of Fig. 9. In 3D, the angle ρ varies between 23.9 • and 24.1 • .
As shown in Fig. 10, Eq. (2) provides an excellent fit of the mean force as a function of time even at the end of the simulation where the tool is close enough to the right wall for the level of the top surface to begin to increase. This happens around 8 s. The fit, according to the Eq. (4), works quite well also in the 3D case, as can be seen in Fig. 11. The parameters values are obtained as in 2D, it yields φ = 18 • , δ = 17 • and B = 4.5 cm. Although the fit in the 3D case slightly underestimates the mean force, the trend is globally well represented and the fitting curve remains within the range of force fluctuations around the mean force. This slight under-estimation can be a finite-size effect in the 3D geometry where there are less particle layers in the simulated system.   We also checked the influence of the cutting velocity on the mean force. At 50 cm/s, we observed a dynamic regime where the granular material intermittently gets loose from the tool, a fact revealed by the force vanishing at some moments. At such high velocities, we observed also a certain fluidization of the free surface manifested by an increase of the void ratio e. Here, we do not enter the details of these observations. Let us only mention that the mean force seems to increase quadratically with the cutting velocity as expected from Bagnold's model [24,25] and found experimentally by Stafford [26].

Force fluctuations
In this section, we will focus on the statistical analysis of variations of the horizontal component of the force supported by the tool around the mean force. We consider two quantities: the deviation ∆F of the force from the mean force and jumps S defined as the difference between two successive extrema of the force. Both of these quantities can be positive or negative. A positive value corresponds to the "loading" of the tool and negative value to its "unloading". We find, however, that the mean fluctuation ∆F and the mean jump S are zero for weak cutting velocities. Moreover, the corresponding histograms are symmetric about the mean value. Figure 12 shows the probability density P (∆F ) of the ∆F for sample A, where we have normalized the fluctuations with respect to the mean value |∆F | of the absolute values of the fluctuations. We do not have enough statistics for the fluctuations beyond ∆F = 7 |∆F | . The largest fluctuations are, however, as wide as 40 |∆F | with the sample D, for example, in which the grains are the biggest. For the force fluctuations not wider than ∆F = 5 |∆F | , we have a large number of events so that the statistical error is too small to be plotted in the figures.
A bias of ∆F and S towards positive values appears as the velocity is increased by one order of magnitude. In this paper, we do not consider such large velocities. Hence, from now on we will no more separate negative and positive fluctuations. Figure 13 shows the probability density P (∆F ) of force fluctuations for the four 2D samples. We see that the four curves collapse on the same curve. We distinguish two different parts on the histogram: (1) For ∆F < 0.5 |∆F | , the histogram is nearly flat; (2) Otherwise, it is a decreasing exponential function: In our simulations, we get α ≈ 1. We remark that, for these statistics, we have removed the end of the simulation, i.e. t > 8 s, since the fluctuations are enhanced as soon as the shear zone touches the right wall. Indeed, we checked that, if the end of the signal is included in the statistics, the tail of the histogram will no more fit to the exponential form (7). In 3D, we find that the histogram is better fitted by a stretched exponential function for large ∆F : with α ≈ 0.4 and β ≈ 1.8. It is remarkable that, while the mean force is the most probable force, the distribution of force fluctuations has the same features as the distribution of contact forces inside a granular medium [12]. Figure 14a,b shows two different states of contact forces in the whole system. Let us notice that a high force acting on the tool corresponds to strongly loaded chains of force near the tool tip, which disappear in case of "unloading" (∆F < 0). The connection between the highly heterogeneous transmission of forces in the contact network, on one hand, and the highly fluctuating behaviour in time, on the other hand, is, however, not straightforward because of force correlations which may extend over distances as large as 10 grain diameters. Other mechanisms can be at the origin of fluctuations. As the Fig. 14a,b. Force networks when the blade is "unloaded" (a) and "loaded" (b) tool cuts the material, the part of the material influenced by the motion of the tool is constantly renewed. Shear lines are thus formed and destroyed. We have not yet investigated this dynamics of shear zones, but we think that it may explain the low frequency fluctuations. Figure 15 shows the distribution P (S) of force jumps S for the four 2D samples. Again, we have an excellent collapse of all data on the same curve when the force jumps S are normalized with respect to the average of the absolute values of S in each sample. In Fig. 15, we also see that a stretched exponential function provides an excellent fit of the data over the whole range of the S: where α ≈ 1.6 and β ≈ 0.8 in 2D. In 3D, we observed the same behaviour with α ≈ 0.5 and β ≈ 1.6. How the force levels ∆F and the force jumps S are correlated? It is obvious that a large negative jump may not occur unless the force level just before the jump is well above the mean force. In the same way, a large positive jump occurs only when the force level just before the jump is far below the mean force. These features can be observed in Fig. 16, where we have plotted the jumps against the force levels from which they start. By definition, |S| should not exceed two times the maximal fluctuation ∆F max in absolute value. Figure 16 shows that most of the ∆F lie within the interval [−5 |∆F | ; 5 |∆F | ] while ∆F max ≈ 20 |∆F | . In this range, in except to a few points, all the S happen to be in the interval [−10 |∆F | ; 10 |∆F | ]. This hints at a stronger correlation than what one could expect from ∆F max ≈ 20 |∆F | . In fact, the lower limit of force jumps as a function of force levels is a well-defined line with slop −1. This is nearly the same for the upper limit. In other words, for a given ∆F , S can have any value in the interval [−5 |∆F | − ∆F ; 5 |∆F | − ∆F ].
A crucial point about the force fluctuations is the way the average values ∆F and S of fluctuations and force jumps vary with grain size. In 2D, where we have four different particle sizes (see Table 1), we expected ∆F to decrease by a factor 1/ √ 3 = 0.58 from sample D to sample A since the mean grain size in A is 3 times smaller than in D. In other words, the number of particles in contact with the tool increases by a factor 3. But we observe only a 10% decrease. As we mentioned above, this should result from space correlations, observed directly on the force map in the form of highly loaded "chains". Of course, these correlations should not be neglected at the particle scale.

Conclusion
In summary, we performed contact dynamics simulations of a noncohesive granular bed containing as many as 10000 particles cut by the motion of a tool. In this paper, we presented data concerning the force exerted by the medium on the tool in the course of its motion. Our analysis of the evolution of the mean force shows that, in spite of large force fluctuations, the data are well fitted by the passive solutions of a limit analysis with Coulomb's failure criterion both in two and three dimensions. This numerical observation thus provides a new confirmation of the relevance of rigid-plastic models for the prediction of the yield stresses of a granular material under the action of a tool.
On the other hand, we found that force fluctuations have a large exponential distribution, independently of the particle sizes. Although it is clear from the simulations that these fluctuations have something to do with the frequent reorganization of the force network and the shear zones, we have yet no quantitative analysis of the microscopic mechanisms at their origin. The presence of large force variations, whose magnitude may considerably exceed the mean force, as computed by a limit analysis, is important for the design of soil-cutting tools. For this reason, a quantitative understanding of the microscopic phenomena in granular materials seems now crucial in many industrial applications. Although the discrete element simulation of a granular assembly of several thousand particles can by no means replace a finite element calculation, the results presented in this paper show that, beside the detailed microscopic data available from such simulations, one may get useful information about the overall behaviour as well, through a suitable statistical analysis of the data.