# Dense granular flow around a penetrating object:

Experiments and hydrodynamic model

###### Abstract

We present in this Letter experimental results on the bidimensional flow field around a cylinder penetrating into dense granular matter together with drag force measurements. A hydrodynamic model based on extended kinetic theory for dense granular flow reproduces well the flow localization close to the cylinder and the corresponding scalings of the drag force, which is found to not depend on velocity, but linearly on the pressure and on the cylinder diameter and weakly on the grain size. Such a regime is found to be valid at a low enough “granular” Reynolds number.

###### pacs:

45.70.-n, 45.50.-jIntroduction.– Describing the motion of an obstacle through granular material is the subject of recent and intensive research with applications to industrial configurations but also to biological and earth science with, example include animal locomotion in sand Maladen2009 () and impact cratering Katsuragi2007 (). If the motion of an object in a simple fluid is known for a long time, especially in the viscous regime at low Reynolds number where the fluid flow and the drag force are analytically known since Stokes’ calculation, the motion of an object in granular matter is still an open question. Such a problem if of fundamental interest, with numerous open questions of statistical physics concerning (for instance) the solid-liquid or jamming transition Dauchot2009 (). Numerous studies have been done concerning the drag force on an object in vertical or horizontal motion in granular matter Albert1999 (); Albert2000 (); Albert2001 (); Stone2004 (); Hill2005 (); Swinney2007 (); Zhou2007 (); Caballero2009 (). All these studies find that the drag force does not depend on the velocity at low velocities, is proportional to the size of the object, and to its depth. As in hydrodynamics, the drag force has been shown to depend on the exact shape of the object Albert2001 (), and also vertical lift forces can develop during horizontal motion Ding2011 (). Flow observations of grains have been also reported in chute flow around a fixed cylinder Kellay2001 (); Zenit2003 () and in the two-dimensional situation of a disk pulled at a constant force in a horizontal assembly of disks on a vibrated plate Dauchot2009 (). Fluctuations have been observed in the force or in the velocity with some “stick-slip” behavior in some cases Albert2000 (); Dauchot2009 (), and the force may depend crucially on the packing volume fraction Swinney2007 (); Dauchot2009 (); Caballero2009 (). In this Letter, we investigate by Particle Image Velocimetry (PIV) measurements the flow around a cylinder penetrating into a dense granular packing together with force measurements. By a continuum hydrodynamic model based on the kinetic theory extended to dense granular systems, we recover the experimental results of the shear localization close to the cylinder with the view of a “hot” cylinder in motion in a viscosity-dependent temperature fluid, together with the good scaling for the drag force.

Experiments.–The experiments consist in a horizontal steel cylinder of diameter 10 mm40 mm and length mm penetrating in a rectangular box (0.1 m in height, 0.2 m in width and 40 mm in thickness) filled with monodisperse millimetric glass beads of diameter 0.5 mm4 mm and density kg m. The granular medium is prepared by gently stirring the grains with a thin rod and the surface is then flatten using a straightedge. We have checked that this preparation leads to reproducible results with only small variations. The solid volume fraction is characteristic of a dense granular packing and the density of the granular medium is thus kg m. The cylinder which is fixed and related to a force transducer by a vertical thin rod, is first above the grain surface and penetrates gradually into the granular packing as the box is raised up by a stepper motor along a vertical translation guide at a constant velocity ranging from 0.1 to 100 mm s. Very careful alignment is taken to prevent any blockage of the cylinder during the motion, and the force at the glass walls without grains is totally negligible. The front and rear wall of the box are in glass allowing visualization of the granular flow around the cylinder. The images taken from a fast video camera (up to 1000 images per second in the full resolution pixels) are analyzed by a PIV software to get the velocity field of the grains. As the camera is fixed in the laboratory frame together with the cylinder, the obtained velocity field shown in Fig. 1 is the velocity of the grains in the frame of reference of the cylinder.

The measured drag force on the cylinder is observed to increase with the depth during its penetration (Fig.2a) with a ratio constant to over the range mm (Fig.2b). The force is found proportional to the cylinder diameter (Fig.2b) and roughly independent of the velocity (Fig.2c). We find also a non-linear dependence of the force with the grain diameter: The force is about constant at large enough grain size ( 1 mm) but increases with decreasing grain size ( 1 mm)(Fig.2d).

Concerning now the grain flow field around the totally immersed cylinder, the first important and key observation is that the velocity field is stationary during the penetration process, meaning that it depends neither on the depth nor on the increasing granular pressure. The velocity field can thus be averaged during each penetration run to extract at each point the mean local velocity whose time fluctuations can be related to the local granular temperature , by . Due to the geometrical configuration, we have used the most appropriate cylindrical coordinates where is the radial distance from the cylinder center and is the angle relative to the downward -axis of motion (with thus down): , with the radial and azimuthal components of the velocity and . As in the classical hydrodynamics situation of a Newtonian fluid, we have checked that and can be decomposed into cosine and sine functions of and radial functions and : and . The radial functions and extracted from measurements in the azimuthal range are displayed in Fig. 3 and show a strong shear localization when compared to the classical viscous Newtonian case, with exponential variations scaled by the cylinder diameter. At a distance from the cylinder surface larger than about one cylinder diameter ( 40 mm in Fig. 3), the grain velocity vanishes. The overshoot of expected from mass conservation in the bidimensional configuration is localized close to the cylinder and relaxes to the asymptotic value 1 (no grain flow) with an inflexion point in the present case in contrast with the long-range decay of the Newtonian case at Re = 0 (Couette-Poiseuille form).

A typical radial profile of the granular temperature along the vertical downward line ( = 0) is displayed in Fig. 4: A domain of roughly constant temperature can be seen around the cylinder with the plateau value followed by an exponential decrease. The temperature profile is roughly the same for different values of and the plateau value is found to vary with the penetration velocity as (see inset of Fig. 4).

Hydrodynamic model.– The granular mean velocity and its fluctuation show important spatial variations. The strong localization of the granular flow that we observe is a rather usual feature of disordered matter where shear bands are commonly observed (see. Schall2010 () for a recent review). We show here that the observed shear bands may be understood using an hydrodynamic description. The starting points are that local momentum and energy balance equations, written as Luding2009 ():

(1) | |||||

(2) |

where is the effective fluid density, the stress tensor, the heat flux, the velocity-gradient tensor, and the temperature-loss coefficient. The stress and the heat flux are related to the velocity, temperature and pressure by phenomenological equations. The choice of these phenomenological equations to describe granular matter is still matter of debate. If momentum transfer and dissipation occur during binary collisions between grains, granular material may be treated using an inelastic gas theory Jenkins1983 (), leading to a Newtonian fluid with , and Fourier’s law for the heat flux. Here is the pressure, the viscosity, the thermal conductivity and the unit tensor. For simplicity, we neglect dissipation and heat transport associated with compressibility and we treat granular matter as an incompressible fluid as experimentally observed. A priori, in a dense granular flow, non-binary collisions cannot be neglected and transport coefficients from inelastic gas theory Jenkins1983 () are no longer valid. By analogy with glassy materials, modification of the viscosity divergence near jamming has been suggested Bocquet2001 (). Numerical simulations of 2D granular materials seem then to show viscosity divergence at packing fractions lower than random close packing Garcia2006 (); Luding2009 (). Since in our experiment, most of the shear is located in a region of high velocity fluctuations, we are not very close to random close packing, and therefore we use the Enskog expressions of the phenomenological coefficients. We obtain Jenkins1983 ():

(3) | |||||

(4) | |||||

(5) | |||||

(6) |

where , , and are non-dimensional functions of the solid volume fraction , and is the velocity restitution coefficient. For , those functions vary with the same dependence on Jenkins1983 (). So we have , and , with , and Jenkins1983 () with a standard value for glass beads. Those expressions of transport and dissipation coefficients emphasize the fact that is not fixed in our experiment but may vary slightly from point to point in response to pressure and granular temperature variations.

The momentum and heat equations (1,2) with dependent transport coefficients are solved numerically for the stationary flow around a cylinder located at the center of a square box. The momentum equation is solved using a Lattice-Bolzmann solver (BGK based D2Q9 model) with non-slip velocity conditions on the cylinder, constant pressure at the top of the square box, and constant upward velocity on the other sides. The heat equation is solved using a finite difference scheme with the condition on the cylinder and on the sides of the square box in agreement with experimental findings. The transport coefficients are taken initially homogeneous in the box and the velocity field is first computed by solving momentum equation with these initial values. From the obtained pressure and velocity fields, the source of heat is calculated and the heat equation is then computed leading to a new temperature field . With the new corresponding fields of transport coefficients as inputs, the momentum equation and then heat equation are solved again and so on. After a few iterations, stationary temperature, pressure and velocity fields which verify (1,2) are obtained. In order to prevent numerical instabilities, viscosity is kept in a finite range. We have checked that the stationary solution is not sensitive to the initial guess of temperature, neither to the cut-off of viscosity. We have also checked that other boundary conditions at the cylinder (partial velocity slip and Robin condition for the temperature) do not change significantly the velocity and temperature fields.

Results.–Our hydrodynamic model reproduces quite well the experimental features observed in the experiments. The radial temperature profile shows the same shape as in the experiments, with a plateau of high temperature close to the cylinder followed by an exponential decrease (Fig. 4). The temperature plateau is also found to be proportional to as in the experiments (see inset of Fig. 4). The computed velocity field is found close to the experimental one as shown in Fig. 3 with a shear localization near the cylinder, and is far from the classical Newtonian case. This shear localization may then be interpreted as a consequence of the strong coupling between viscosity and temperature. At a given pressure the viscosity varies as . So, at a given viscous shear stress , the production of heat is proportional to . This mechanism creates a self-lubricating layer of low viscosity near the cylinder.

The drag force on the cylinder is calculated by integrated the stresses on the disk, taking into account both the “pressure” term (from normal stresses) and “viscosity” term (from shear stresses). If these two terms are equal in the Newtonian case, the pressure term is here about twice the viscosity term. The total calculated drag force is found to be independent of the velocity, proportional to the cylinder diameter and to the pressure, in agreement with the experimental observations by considering that pressure is proportional to depth. That dependance may be understood in the hydrodynamical model by considering a “granular” Reynolds number Re = based on the viscosity near the cylinder. We have checked that all these previous findings correspond to a low Reynolds number regime (Re ). In low Re hydrodynamics, one expects that the force scales as , and then varies here linearly with the pressure, and independently of the velocity as is proportional to the pressure and to thus to . We also found numerically the same non linear variation of the drag force with the grain size as in the experiments. Note that this variation is hard to infer simply from the set of equations. When Re , the velocity field no longer exhibits up/downstream symmetry, and the pressure and temperature profiles are also different from the Re case.

Concluding remarks.–We have investigated experimentally the penetration of a cylinder at a constant velocity inside a dense granular packing with both force and velocity field measurements, and we have modeled this problem by a continuum hydrodynamic approach. The finding of a strong shear localization close to the cylinder can be viewed by the coupling between viscosity and temperature in the problem of a self-heated cylinder. The localization of the flow near a sedimenting hot sphere has indeed already been reported in classical fluids with temperature dependant viscosity Ansari1985 (). Such a shear localization has been seen also for a sphere sedimenting in a non Newtonian fluid with shear thinning behavior Atapattu1995 (). The experimental findings of a force regime independent of velocity, proportional to the depth and to the cylinder diameter are recovered by our model based on kinetic theory adapted for dense granular systems, and have been shown to correspond to an hydrodynamic regime of low “granular” Reynolds number. Other models exist for dense granular flows such as the one based on a local rheology with a friction coefficient depending on a non-dimensional shear rate Pouliquen2006 (). Such models may also lead to the observed flow localisation since it corresponds to visco-plastic/shear thinning behavior. Another issue to explore would be the different force value measured in the plunging and the withdrawal situations Hill2005 (); Swinney2007 (), with the role played by gravity and the boundary conditions (bottom wall vs free surface). In addition, the force scaling is expected to change in a higher Reynolds regime from a “viscous” scaling to an “inertial” scaling which may explain the complex force terms measured in impact situations Katsuragi2007 (); Goldman2008 (); Seguin2009 () where the “granular hydrodynamic” regime change certainly from high to low Re during the penetration process. Applications to non-stationary granular flows may also been considered.

We thank J.T. Jenkins, E. Trizac and E.J. Hinch for stimulating discussions. We are grateful to F. Martinez for its contribution to the experiments and A. Aubertin for the development of the experimental setup. This work is supported by the ANR project STABINGRAM No. 2010-BLAN-0927-01.

## References

- (1) R. D. Maladen et al., Science 325, 315 (2009).
- (2) H. Katsuragi and D. J. Durian, Nature Physics 3, 420 (2007).
- (3) R. Candelier and O. Dauchot, Phys. Rev. Lett. 103 128001 (2009).
- (4) M. Schröter et al., Europhys. Lett. 78, 4404 (2007).
- (5) G. A. Caballero-Robledo and E. Clément, Eur. Phys. J. E 30 395 (2009).
- (6) I. Albert et al., Phys. Rev. Lett. 84, 5122 (2000).
- (7) R. Albert et al., Phys. Rev. Lett. 82, 205 (1999).
- (8) I. Albert et al., Phys. Rev. E. 64, 061303 (2001).
- (9) M. B. Stone et al., Nature 427, 503 (2004).
- (10) G. Hill, S. Yeung, and A. Koehler, Europhys. Lett. 72, 137 (2005).
- (11) F. Zhou, S. G. Advani, and E. R. Wetzel, Phys. Fluids 19, 013301 (2007).
- (12) Y. Ding, N. Gravish and D. I. Goldman, Phys. Rev. Lett. 106, 028001 (2011).
- (13) Y. Amarouchene, J. F. Boudet and H. Kellay, Phys. Rev. Lett. 86, 4286 (2001).
- (14) D. Chehata, R. Zenit, and C. R. Wasgreen, Phys. Fluids 15, 1622 (2003).
- (15) P. Schall and M. van Hecke, Annu. Rev. Fluid Mech. 42, 67 (2010).
- (16) S. Luding, Nonlinearity 22, R101 (2009).
- (17) J.T. Jenkins and S.B. Savage, J. Fluid Mech. 130, 187 (1983).
- (18) L. Bocquet et al., Phys. Rev. E 65, 011307 (2001).
- (19) R. García-Rojo, S. Luding and J.J.Brey, Phys. Rev. E 74, 061305 (2006).
- (20) A. Ansari and S. Morris, J. Fluid Mech. 159, 459 (1985).
- (21) D.D. Atapattu, R.P. Chhabra, and P.H.T. Uhlherr, J. Non Newtonian Fluid Mech. 59, 245 (1995).
- (22) P. Jop, Y. Forterre, and O. Pouliquen, Nature 441, 727 (2006).
- (23) D. I. Goldman and P. Umbanhowar, Phys. Rev. E 77, 021308 (2008).
- (24) A. Seguin et al., Europhys. Lett. 88, 44002 (2009).