Dense flow around a sphere moving into a cloud of grains

A bidimensional simulation of a sphere moving at constant velocity into a cloud of smaller spherical grains without gravity is presented with a non-smooth contact dynamics method. A dense granular “cluster” zone of about constant solid fraction builds progressively around the moving sphere until a stationary regime appears with a constant upstream cluster size that increases with the initial solid fraction φ0 of the cloud. A detailed analysis of the local strain rate and local stress fields inside the cluster reveals that, despite different spatial variations of strain and stresses, the local friction coefficient μ appears to depend only on the local inertial number I as well as the local solid fraction φ, which means that a local rheology does exist in the present non parallel flow. The key point is that the spatial variations of I inside the cluster does not depend on the sphere velocity and explore only a small range between about 10−2 and 10−1. The influence of sidewalls is then investigated on the flow and the forces.


Introduction
The response of a granular material to a mechanical perturbation arising from the motion of a solid object at its surface or in its bulk is a fundamental issue in many fields such as civil engineering for the resistance of soils to the penetration or extraction of stakes and piles, biophysics for the understanding of animal locomotion in sand which may inspire new robotics, geophysics for the collision phenomena in the impacts of meteorits on planets and asteroids, or in the formation of protoplanetary disks from dust particles.In all these cases, the complex rheology of the granular material plays a key role and needs to be wellunderstood, which is particularly hard when the packing is dense with a high solid fraction φ close to the so-called liquid/solid or jamming transition.Different rheological laws have been proposed for dense granular flows, such as the so-called μ(I) rheology which relates the local friction coefficient μ and the local packing fraction φ to the dimensionless inertial number I [1].Such a rheological framework gives good results in quasi-steady and quasiparallel flows but its validity for strongly non parallel flow still needs to be addressed.Even in simple parallel shear flows, the existence of such a local rheology and the influence of solid walls is still under investigation [2].As a matter of fact, the presence of far boundaries such as fixed or mobile solid walls has shown to have strong and non elucidated actions in different situations [3,4].To examine this important question, one must have access simultaneously to the strain and stress fields in the bulk flow which e-mail: antoine.seguin@u-psud.fr is possible experimentally in some cases [5] but now easier with discrete numerical techniques.
The object of the present paper is to investigate the local rheology of granular matter in a strongly non parallel flow around a moving sphere without any gravity field.We use for this a recent powerful numerical method to investigate such a sphere motion within a bidimensional cloud of spherical particles initially at rest with a solid fraction φ 0 far below the jamming point.We first recall the main results reported in [6] in the case where the flow occurs far from any boundaries, before discussing the possible influence of sidewalls.

Flow configuration and numerics
The numerical method we use here belongs to the class of "Non-Smooth Contact Dynamic methods" [7] and is described in details in [8].We simulate the motion of dissipative rigid spheres with a non-elastic impact law (zero restitution coefficient for the collisions) but without any static nor dynamic friction between the spheres.The granular assembly, made of slightly polydisperse spherical grains of mean diameter d g and density ρ to avoid any possible crystallization, is contained in a rectangular box of size L x × L y with (x, y) = (0, 0) at the center of the box.The aspect ratio of the box has been varied in the range 1 L x /L y 15.An intruder sphere of larger diameter d = 10d g , which is initially placed at one side of the rectangular box at x −L x /2 and y = 0, is then moved at constant velocity in the x-direction toward the opposite side.The confinement by the sidewalls has been varied The interaction between the grains and the moving sphere as well as between the grains and the walls is solved using the same law as for the grain-grain interaction.To prepare the granular medium at the initial solid fraction φ 0 = N(πd 2 g /4)/(L x L y − πd 2 /4), we pick randomly N centers of non-intersecting spheres of mean diameter d g and apply a random force on the grains to ensure an uniform spatial configuration of spheres without any contact.All the explored initial configurations are far below the jamming point φ J 0.84.
Starting the intruder sphere motion, the numerical solver gives the velocity v m of each particle m and the contact force f nm exerted by particle n on particle m, which allows to calculate the corresponding stress tensor σ m = (4/πd 2 g ) n e mn ⊗ f nm , where e mn is the unit vector giving the direction from the center of particle m toward the center of particle n and ⊗ is the vector outer product.The local solid fraction φ m around each particle m is computed using a Voronoi tesselation.The local solid fraction φ at any point of the domain is then computed by interpolating the φ m values onto a cartesian grid of step Δx = Δy = 0.05d g , as well as the local velocity field v and stress field σ from the v m and σ m values.From the rate of strain tensor From the stress tensor, we extract the pressure p = −(1/2) k σ kk and the shear stress τ = [(1/2) i, j (σ i j + pδ i j ) 2 ] 1/2 .The drag force F exerted by the grains on the intruder can be calculated as F = − m f m0 • e x , where f m0 is the contact force exerted by particle m on the intruder and e x is the unit vector along the x direction of motion.
The spatial variation of the different quantities will be expressed in the polar coordinates (r, θ), where θ is the angle relative to the x direction of motion and r the radial position relative to the moving sphere center.

Flow results without sidewalls
As the intruder starts moving, it perturbs progressively the initial grain assembly and a dense cluster zone of touching  grains grows around as illustrated in Fig. 1.The position r c (θ, t) of the corresponding front delimiting the inner perturbed cluster zone to the outer unperturbed zone is extracted from the simple criterion p(r > r c ) = 0 for each θ.A triangular zone with no grain inside, thus appearing in white in Fig. 1, exists in the wake of the intruder as already reported and analyzed in details by [9] in bidimensional experiments.In the present paper, we only focus on the cluster zone upstream the intruder (−π/2 θ π/2) which is the key region where the drag force originates, with high stresses and strain rates inside.The upstream extension of the cluster λ = r c − d/2 is averaged in the small θ range −5 deg < θ < 5 deg around the x-direction of motion.As in one dimensional experiments where a straight rake starts moving [10], we observe that the front position λ moves away linearly in time, with a velocity that is proportional to V 0 and increases with φ 0 .But in contrast to the one-dimensional configuration of [10], this regime is here only transient until a steady regime is reached with a constant value λ c despite some fluctuations (Fig. 2).This steady regime appears when the grain flux upstream the intruder is balanced by the grain flux on the intruder sides which can not occur in the one-dimensional configuration of [10] as no grains can circumvent the rake.In the following, we restrict our study to this steady state regime which allows time averaging of the different quantities.
The cluster steady size λ c is observed to increase with φ 0 from only a few grains at low φ 0 (e.g.λ c 3d g at φ 0 = 0.3) to many grains at high φ 0 (e.g, λ c 53d g at φ 0 = 0.7), and diverges at the approach of a critical value φ c with the scaling law λ c 1.5d g (φ c − φ 0 ) −2 with the value φ c 0.85 ± 0.01 very close to the jamming point φ J .The local solid fraction φ in the cluster does not vary significantly with θ in a large azimutal range and displays the about constant value φ 0.83 in all the cluster for large enough φ 0 (φ 0 0.6), except close to the intruder where φ decreases down to about 0.75 and near the front where φ decreases down to φ 0 .
In the present simulations where no gravity acts and no pressure is imposed from any external boundary, no stress scale exists except the kinetic pressure ρV from collision processes.The present regime corresponds therefore to the inertial high velocity regime found by [11] in their bidimensional experiments of a disk dragged within a monolayer of steel beads and also found by [12] in their numerical simulations of a bidimensional assembly of disks with no friction with the bottom plate.In our simulation, we do not observe any quasi static regime at low velocity where the drag force would be velocity independent or would depend only weakly on the velocity.The existence of such a quasi-static regime which is the most often seen regime [3,4] arises from the existence of another natural scale of pressure in the system which may come either from gravity and solid friction.We found here that the drag force F increases linearly with the initial solid fraction φ 0 with the scaling law F 5φ 0 ρd g 2 V 2 0 φ 0 ρdd g V 2 0 /2.Thus, F does not diverge at the approach of φ c in contrast to the cluster size λ c .This may be explained by the fact that the grains at the outer limit of the growing cluster do not touch any limit boundary in the present simulations.
Let us now detail the flow in the dense cluster around the moving intruder (Fig. 3).The dilation rate ε is about zero and thus the flow is incompressible in the overall cluster, except near the cluster rim at the front (Fig. 3a) where there is significant φ variation collinear with velocity as mass conservation reads φ ε + v • ∇φ = 0.The shear rate γ(r, θ) is observed to be maximum close to the moving sphere at r d/2 and then strongly decreases away from it with the law γ(r, 0) 4(V 0 /d g )[(d g /r) 2 + α] with α 0.003 (Fig. 3a), and also decreases toward the equator.The key point is that its maximal value γ(d/2, θ) does not depend significantly on φ 0 neither does its radial or azimuthal decreasing rate.This means that an intrinsic flow appears close to the moving sphere within the cluster zone independently of its possible diverging size.This flow has an intrinsic length scale which is thus independent of φ 0 and is of the order of 1d and thus of about 10d g for the present size ratio d/d g = 10.This flow zone corresponds  to the narrow crown of lower solid fraction φ < 0.8 (Fig. 3c).
The local pressure p(r, θ) is maximal at the sphere surface in front of the moving sphere (r d/2, θ 0) and decreases radially away from it with the approximate law p(r, 0) 5.5φ 0 ρV 2 0 [(d g /r) + β] with β 0.07 and toward the equator (Fig. 3b).The spatial variations of the shear stress τ are very similar to those of the pressure p with always τ < p (Fig. 3b).
The strong coupling between τ and p means that the rheological behavior of the grain assembly appears to be of frictional type although no microscopic friction exists in the system.We thus test the possible existence of the local rheology μ(I) where the local friction coefficient μ = τ/p would be linked to the local inertial number I = γd g / p/ρφ [1].As τ and p have the same scaling in V 0 , μ will not depend here on the velocity V 0 neither does I values as both γ and √ p scales as V 0 .The only spatial variation of I will come from the weak difference of the spatial scalings of γ and √ p, together with the weak spatial variation of φ (see φ(r) shown in Fig. 3d).We see in Fig. 3d that I has its highest value of about 0.2 close to the intruder and decreases away from it following the same master curve whatever the cluster size.For the largest cluster size (φ 0 = 0.75), we see that I tends towards a non-zero asymptotic value close to 10 −2 .We observe only weak azimuthal variations of I around the intruder.In Fig. 4 where μ and φ are plotted as a function of I, all data collapse into master curves which means that a local rheology appears in the present flow.The solid fraction φ decreases with the approximate scaling φ = φ m (1 − aI) with φ m = 0.846 ± 0.002 and a = 0.80 ± 0.03, and the friction coefficient μ increases as μ = μ s + (μ 2 − μ s )/(1 + I/I 0 ) with μ s = 0.25 ± 0.05, μ 2 = 0.78 ± 0.02 and I 0 = 0.03.The φ m value is very close to the critical value φ c 0.85 found for the divergence of λ c .We believe that φ c and φ m both correspond to the jamming point φ J of the system.The about constant value 0.83 observed for φ in the cluster core corresponds to the asymptotic minimum value of I m = α/β 1/2 0.01 reached in the flow.The observed μ(I) In the present configuration, the inertial number I that naturally arises ranges within only one decade from 0.01 to 0.15.As a consequence μ only varies roughly from 0.45 to 0.70 and φ roughly from 0.75 to 0.83.The region of high I (high μ and low φ) is close to the moving sphere but it is worth noting that I does not vanish far away.Note that the present natural range of I corresponds to a well-posed behaviour of the μ(I)-rheology [14].

Influence of sidewalls
Let us now present some preliminary results on the influence of sidewalls.By decreasing transverse wall gap L y from 800d g to 40d g , we observe that the drag force F increases with the confinement ratio d/L y as expected from [13], and that the upstream cluster size λ c also increases.When the sidewalls are close enough(L y /d 10, the front of the cluster is no more curved but perpendicular to the sidewalls and much far from the intruder.If we look at the local shear rate and pressure (Fig. 5), we observe that γ is not affected significantly upstream the intruder by the presence of these lateral walls but the pressure is, with a non-zero plateau value that appears far from the intruder in the cluster beyond the rapid pressure decrease close to the intruder.An approximate scaling for the pressure is now p(r, 0) 5.5φ 0 ρV 0 2 [d g /r + β(d/L y )], where β(d/L y ) is a increasing fonction of d/L y that still needs to be determined.

Conclusion
Our simulation results show that a steady state regime is reached in the dense flow around an intruder moving at constant velocity within a cloud of grains.The granular flow in the dense cluster around the intruder is found to obey a local rheology where the local friction coefficient μ and the local solid fraction φ are given by the inertial number I even if no microscopic friction between the grains is considered here.With close sidewalls, the spatial variation of stresses are strongly affected with a significative non-zero value that arises for the pressure far upstream the intruder within a long range plateau.

Figure 1 .
Figure 1.Snapshot of a simulation for the initial solid fraction φ 0 = 0.75.The grains appear in (light) dark blue when there is (no) contact force in between.The larger intruder moving from left to right appear in black.The length λ c corresponds to the cluster size in front it.

Figure 3 .
Figure 3. Radial evolution of (a) the dimensionless dilation rate ˙ and shear rate γ, (b) the dimensionless pressure p and shear stress τ, (c) the solid fraction φ, and (d) the inertial number I for an initial solid fraction φ 0 = 0.75.