## Abstract

Experiments are conducted to measure the resistance experienced by light cylinders rolling over flat beds of granular media. Sand and glass spheres are used for the beds. The trajectories of the rolling cylinders are determined through optical tracking, and velocity and acceleration data are inferred through fits to these trajectories. The rolling resistance is dominated by a velocity-independent component, but a velocity-dependent drag exceeding the expected strength of air drag is also observed. The results are compared to a theoretical model based on a cohesionless Mohr–Coulomb rheology for a granular medium in the presence of gravity. The model idealizes the flow pattern underneath the rolling cylinder as a plastically deforming zone in front of a rigidly rotating plug attached to the cylinder, as proposed previously for cylinders rolling on perfectly cohesive plastic media. The leading-order, rate-independent rolling resistance observed experimentally is well reproduced by the model predictions.

## 1. Introduction

It is commonly experienced that motion over loose sand is hindered by the response of the granular medium to the applied load. In particular, sandy surfaces often present substantial difficulties to wheeled vehicles [1], impacting the design of tyres and robotic rovers [2,3]. The rolling of rocks over sand may also play an important role in geophysical processes such as the formation of screes [4]. Although Coulomb provided detailed measurements of rolling resistance, Reynolds [5] gave the first quantitative discussion of the origin of the phenomenon in terms of the deformation of the surface in rolling contact. The resulting friction experienced by rubber and metal surfaces is now understood to be a combination of Reynolds’ deformation and surface adhesion [6,7]. For wheels travelling over sand, the problem is complicated further by the lack of a general theory of the deformation of a granular medium, the significant depth to which the wheels of a heavy vehicle will sink and the substantial sideways flow of sand. Nevertheless, a rich vein of literature has built up on the resistance to wheels driven over sand and soft soils (e.g. [1,3,8–16]).

In this work, we consider the rolling motion of a long and light cylinder over an initially static granular bed. In this limit, the contact area between the cylinder and the bed is small relative to the surface area of the cylinder, and the deformation of the granular media is largely two dimensional and both localized and superficial. This sets our work apart from much of the literature on rotating wheels and recent experiments on rolling spheres [17], where the problem is inherently three dimensional and the weight of the object traversing the surface leads to significant sinking and compaction of the underlying granular material. Additional complications are thereby minimized, which may allow us to gain better insight into the dynamics of granular deformation under a rolling object.

Our study begins with a series of experiments to determine the resistance experienced by the rolling cylinders. A particular focus is to experimentally examine the dependence of the resistance on the speed of the cylinders, which reflects the rate-dependent rheology of the granular material and bears on the relevance of recently proposed empirical friction laws [18] For this task, we use two types of granular media: a coarse sand and glass spheres (ballotini). For both, the grain size is much smaller than the radii of the rolling cylinders.

We complement the experiments with a rate-independent theoretical model for the leading-order resistance based on Mohr–Coulomb plasticity theory for the granular bed. In the limit that the contact area is much smaller than the cylinder area, the stress field induced by the rolling cylinder can be constructed by asymptotic theory. This device, which was proposed by Spencer [19,20], has been used previously to model the rolling of cylinders over cohesive material [21,22], neglecting the effect of gravity. For a non-cohesive granular medium with a free surface, however, gravity cannot be neglected, as in the calculation of the critical load supported by a foundation [23–25]. Our analysis applies Spencer’s method to a cylinder rolling over such materials.

A critical detail of the generalization is the flow pattern that must be adopted underneath the contact region between the sand and the cylinder. We adopt a construct proposed by Collins [26] (see also [27]) in which the deformation due to rolling consists of two parts: a zone of plastic deformation extending below and ahead of the forward section of the contact region, together with a plug that is rigidly attached to the cylinder below the rear section of the contact region. Collins’ construction allows for contacts that are not necessarily small and has received some support from numerical computations [15], but was proposed for cohesive material neglecting gravity. The two-part structure of the flow pattern is, however, essential for a consistent description of the forces and torques acting on the cylinder, and is reminiscent of experiments and models of wheels rolling over sand [9–14].

The experimental apparatus, procedures and dataset are described in §2. In §3, we outline the theoretical model to calculate rolling resistance. We compare the predictions of the model to our experimental results in §4. The appendices contain some additional theoretical results, exposing issues with the Mohr–Coulomb predictions for the velocity field and free surface of the granular medium, and relating our analysis back to that of Collins [26].

## 2. Experiments

### (a) Apparatus and protocol

We constructed a channel using a high-density polyethylene board and wooden side rails to contain a bed of granular material, over which we rolled a variety of cylinders (figure 1). The channel was approximately 244 cm in length, 31 cm in width and 9 cm in depth. Two granular materials were used in our experiments: PetCo aquarium sand and Potters Industries A-100 glass beads (ballotini). Both materials had mean grain sizes near 1 mm, but the sand grains were more angular than the nearly spherical ballotini, resulting in substantially different dynamic angles of repose between the materials. The material parameters are provided in table 1, and were taken from Sauret *et al.* [28], who used the same materials.

To prepare the granular bed, we filled the channel with a granular material, and then loosely scraped off the top surface layer down to a depth set using the top of the side rails. No effort was made to compact the bed any further. To reset the bed after a disturbance to the surface, we manually stirred the material throughout its depth, and repeated the levelling procedure. On the ballotini, the cylinders typically left large (several millimetres deep) tracks as they rolled, so the surface was reset between each run. On the sand, the cylinders left tracks that were barely perceptible to the naked eye, but could be highlighted by shadows cast when the surface was illuminated from a shallow angle. The shadows revealed that the surface traversed by the cylinders had been smoothed out over the grain scale, and small levees had been left to either side of the track. To speed up the data taking process, we therefore only reset the sand surface when a large disturbance occurred (e.g. a stray cylinder). Despite the lack of any substantial surface rearrangement by the cylinders on the sand, the rolling resistance was measurably sensitive to the surface history: the resistance was noticeably different when rolling a cylinder over previous tracks (an effect familiar from the rolling of spheres over rubber and metal [6]). To reduce this effect as much as possible, we began the experiments with sand using a bedding-in procedure described in §2c.

The cylinders were sections cut from polyvinyl chloride (PVC) and other plastic tubing and had the radii *R* and masses *M* listed in table 2. This table also reports the quantity *α*=*I*/*MR*^{2}, where *I* is the moment of inertia, which provides a dimensionless description of the radial mass distribution, with *α*=1/2 for a uniform solid cylinder and *m* of the bed, providing a resolution of roughly 1 *mm* per pixel.

The cylinders were placed on a flat inclined ramp at the top of the runway and rolled onto the surface of the granular bed. A flexible sheet was placed over the end of the ramp to smooth the transition onto the bed and reduce any bouncing of the cylinders. The initial velocity of the cylinders as they entered the bed was controlled by adjusting the starting position of the cylinders along the ramp. Trials were discarded if they came within a few centimetres of the channel side walls, which tended to occur more frequently for lighter cylinders and at higher initial velocities. We repeated rolls until we recorded three to five successful trials with the same experimental parameters.

### (b) Data processing

The data processing consisted of three stages. First, ‘image trajectories’ indicating the image coordinates of the cylinders for each video frame were calculated (figure 1). Second, position markers on the channel were used to infer the camera’s projection for the scene. Third, the ‘object trajectories’ indicating the real-space positions of the cylinders over time were extracted by inverting the camera projection for the calculated image trajectories. For the latter, we refer to the tops of the side rails to define a laboratory coordinate system, with *X* pointing along the channel, *Y* across it and *Z* directed vertically upwards.

#### (i) Image trajectories

The bright orange paint of the cylinders allowed them to be precisely differentiated from the bed and the rest of the scene via colour filtering. The OpenCV 3.0 library [29] was used to extract individual frames from the videos and convert them to hue-saturation-value (HSV) colour space. We applied a light (11 pixel) Gaussian blur to each frame and computed the HSV colour distance *d* of each of the pixels within the granular bed from a reference orange value of (4 205 255). Using the colour distance *d*, we constructed a ‘match value’ *m* for each pixel as

#### (ii) Camera calibration

To accurately determine the cylinder’s position on the bed from its image coordinates, the camera projection was precisely characterized. We performed a standard calibration using routines from the OpenCV 3.0 library, which are based on a pinhole camera model with corrections for common modes of lens distortion. We first used images of a standard chessboard pattern to determine the camera’s intrinsic projection parameters. For each scene, we then used the image positions of distance markers on the channel side rails to infer the pose of the channel relative to the camera.

The image coordinates of any point (*X*,*Y*,*Z*) in the laboratory frame could then be estimated by composing the best-fit pose transformation with the intrinsic camera projection. The difference between the measured image coordinates of the distance markers and the estimated image coordinates computed via this procedure constitutes the reprojection error. Our calibration procedure typically yielded reprojection errors of approximately 3 pixels.

#### (iii) Object trajectories

The camera and scene calibrations were used to convert the image trajectories into object trajectories describing the (*X*,*Y*) position of the cylinder centres over time for each trial. In a given image, each pixel is the projection of an entire line-of-sight in the object coordinate space. If the image coordinates and one of the object coordinates (e.g. *Z*) of a point are known, however, this degeneracy can be broken and the projection inverted to determine the two remaining unknown object coordinates (e.g. (*X*,*Y*)). For each trial, we used the vertical offset of the material surface from the side rails and the radius of the cylinder to estimate the height *Z* of the plane containing the cylinder centre. The inversion of a cylinder’s image position at this *Z* then yields an estimate of the (*X*,*Y*) coordinates of the centre of the cylinder over the bed. The error in this procedure is dominated by the reprojection error of the camera calibration, resulting in a systematic uncertainty of approximately 3 mm in the absolute positions in our object trajectories.

### (c) Results

We initially tilted the channel by supporting it at different heights under the two ends, and varied the angle to look for steady rolling states. We found, however, that the channel was insufficiently stiff and would flex under its own weight, resulting in a spatially varying downslope component of gravity, as well as further errors in the calibration procedure. To avoid these issues, we laid the channel level with distributed supports to keep it flat to within 1 mm. Removing any gravitational acceleration along the bed had the further advantage of ensuring that the rolling resistance was more easily measured. However, because our results pertain to unsteady motion, they may differ from those for cylinders towed at constant speed.

In our preliminary experiments on the sand, we also explored the effect of the length of the cylinders and the depth of the granular bed on the resistance. We found that there were no significant differences in the trajectories of different length cylinders of the same PVC piping. Similarly, no significant or systematic differences were found for trials on 1, 3, 5, 7 and 9 cm deep sand beds. We therefore decided to fix the cylinder lengths to be 15 cm and the depth of the bed to be 7 cm for the sand and 2 cm for ballotini to reduce the number of parameters being varied. With these considerations, we built a dataset of over 400 trials. Plots of the trajectories for several cylinders on the ballotini are shown in figure 3.

Two-dimensional *X*−*Y* plots of the trajectories are useful for identifying trials in which cylinders swerved to the side and rolled close to the side rails, where they tended to slow down or abruptly change direction. Many trails which strayed from the centre line of the bed appeared to continue turning as they rolled (figure 3*a*), suggesting that there may be an instability associated with perturbations away from straight or level rolling. The veering rate correlated with mean cylinder density: overall, 15% of the trials veered more than 4 cm from the centre line, but 25% of the lighter cylinders and under 10% of the heavier cylinders rolled off centre. To minimize any effects of the side walls on the bed dynamics, we deleted the trials that veered more than 4 cm off-centre, and used the *X* components of the remaining trajectories to study the kinematics of the cylinders.

The millimetre-scale accuracy of the position determination allows us to use finite differences in time to estimate the downslope velocities. Further differentiation to determine the acceleration, however, is precluded by the noise in the position measurements. In figure 3, the cylinder accelerations are instead estimated by differentiating linear least-squares fits of quadratic polynomials to the velocity data in figure 3*c*). In general, however, we find that the form of such fits often severely biases the shape of the resulting acceleration–velocity curves (figure 3*d*), and ultimately, we favour the forward-modelling approach of §4 for determining the dependence of acceleration on velocity.

For trials on the sand, where the surface was not reset between each run, we observed a systematic change in the stopping distances of the cylinders when they were first rolled over a fresh surface. We attributed this to a superficial compaction and smoothing of the top layers of grains, which affected the rolling resistance. To eliminate this effect and begin with a surface that accommodated any sideways migration of our cylinders, we therefore followed a bedding-in procedure at the beginning of the experiments: we first rolled a much longer (25 cm) cylinder down the track repeatedly until it achieved a consistent stopping distance. Despite this preliminary preparation of the surface, there were still noticeable changes in the rolling resistance when cylinders were switched. The stopping distances of the first few trials with each cylinder typically increased or decreased, depending on the cylinder, until saturating approximately 10 cm from the first stopping position (figure 4*a*). After three or more trials, the stopping distances were typically reproducible to within several centimetres, leading us to discard the initial trials showing obvious re-bedding behaviour. Evidently, each cylinder affects the granular surface slightly differently, but we did not explore this effect in any further detail.

Another interesting feature of the rolling dynamics arises at the ends of the trajectories where the cylinders come to a stop on the bed. As shown in figure 4*b*, the cylinders tend to roll backwards several millimetres before coming to a complete stop. This effect may be due to the mound of grains, or ‘bow wave’, that builds up ahead of the cylinder and transmits the forces from the bed during rolling. Once the cylinder is arrested, however, the horizontal bed reaction becomes unbalanced, momentarily pushing the cylinder backwards.

To test the reproducibility of the dynamics, we conducted trials with a range of initial velocities for each cylinder on the sand. For initial velocities in the range of 50–100 cm s^{−1}, the deceleration of the cylinders did not show a pronounced dependence on the initial velocity, as can be seen in figure 4*c*. Trials with higher initial velocities around 175 cm s^{−1}, however, showed substantially larger decelerations. It was difficult to produce consistent trials at these speeds, as the cylinders would frequently bounce after leaving the ramp and veer abruptly into the side rails. We therefore exclude these higher-velocity runs from the bulk of our analysis, but note that the discrepancies in figure 4*c* suggest that efforts to characterize rolling resistance only in terms of rolling speed may well conceal the full nature of the dynamics at high speeds.

Besides the resistance encountered due to the deformation of the granular surface, air resistance may play an important role in some trials. Taking 50 cm s^{−1} as a characteristic velocity, 6 cm as a characteristic cylinder diameter, and *ν*≈1.5×10^{−5} *m*^{2} *s*^{−1} for the kinematic viscosity of air, we estimate the Reynolds number of the airflow around the rolling cylinders to be of order 2000. Thus, air drag on the cylinders is expected to depend quadratically on velocity with an order-1 drag coefficient. To directly measure the effect of air drag, we performed additional trials in which we rolled the cylinders over a glass sheet placed on the channel (see §4).

Finally, several trials were examined to determine whether any sliding occurred as the cylinders rolled over the granular surface. A thin black line segment was marked on the outside of each cylinder, and the times when these lines crossed the top of the cylinders were determined from the videos. The corresponding object positions were then used to measure the distance travelled during one complete rotation. In the cases examined, this distance was indistinguishable from the cylinder circumference to within the experimental precision of approximately 1%, based on our ability to determine the time at which each rotation was completed. Thus, we did not observe any noticeable slip (true or effective; see §4b(iv)) for any trials.

## 3. Modelling

### (a) Cylinder dynamics

To model the dynamics of the system, we consider the motion of a two-dimensional rigid cylinder rolling over a deformable granular surface. The geometry is sketched in figure 5. We define a new cartesian coordinate system centred under the leading contact point at the neutral surface of the sand, with *x* pointing forwards and *z* upwards. The cylinder has radius *R*, mass per unit length *M* and moment of inertia per unit length *I*=*αMR*^{2}. The cylinder has velocity *Ω*(*t*). The translational and rotational dynamics of the cylinder are given by
*g* is gravity, and (*F*_{x},*F*_{z}) and *T* are the contact force and anticlockwise torque from the underlying medium, respectively, again per unit length of the cylinder.

The contact forces are computed by integrating the surface stresses over the contact arc ** σ** is the granular stress tensor, and

**n**and ℓ denote the cylinder normal and arc length along

**r**is the position vector from the bottom of the cylinder.

The horizontal and angular momentum equations can be combined to give

Below, we evaluate the integrals in (3.4) and (3.7) to leading order in the small parameter *a*/*R*, where *a* is the horizontal length of the contact arc, treating the sand as a plastic material satisfying the Mohr–Coulomb constitutive law. Before accomplishing this feat, we note that if the surface stresses are *O*(*ρ*_{s}*ga*), then *F*_{x} and *F*_{z} are formally *O*(*ρ*_{s}*ga*^{2}) from (3.4), while *F*_{r} is *O*(*ρ*_{s}*ga*^{3}/*R*) from (3.7). However, if *α*)*F*_{x}=−*F*_{r}, then *F*_{x}∼*F*_{r}. This can only be arranged if *F*_{x} turns out to contain a small numerical factor of order *a*/*R*, in which case *F*_{x}∼(*a*/*R*)*F*_{z}. For the contact area to remain small during rolling, *O*(*a*/*R*) in comparison to *F*_{z}≈*Mg* in (3.2). This allows us to estimate the magnitude of the horizontal accelerations as *F*_{z}≈*Mg*∼*ρ*_{s}*ga*^{2} and *M*∼*ρ*_{c}*R*^{2}, and hence the contact area scales as

### (b) Sand dynamics

#### (i) Slip-line theory

To form a tractable model, we neglect inertial forces in the bed, which scale like *ϕ* is the internal angle of friction, *p* is the mean compressive stress and *θ* is the anticlockwise angle from the horizontal to the direction of the least compressive stress. With this rheology the stress field is conveniently represented in terms of its characteristic curves, the slip lines. These curves are composed of two families; the *α*-lines are defined by
*β*-lines satisfy

#### (ii) Boundary conditions

We calculate the leading-order forces in the small parameter *a*/*R*, for which it is sufficient to linearize the surface conditions about *z*=0. In the coordinate system centred at the leading edge, the section −*a*<*x*<0 is the contact surface *x*>0 is a free surface.

We impose a friction condition beneath the cylinder that translates to demanding that
*θ*_{f}=0 so that *σ*_{xz}=0 and *θ*_{f}=−*ε* and the *β*-lines are parallel to the surface. If *δ* denotes the surface friction angle, then

The free surface ahead of the contact area is in a state of compression with
*σ*_{xz}=*σ*_{zz}=0. This leads to a subsurface triangular region (a ‘Rankine zone’) where the stress is hydrostatic (*σ*_{zz}=*ρ*_{s}*gz* and *ε* to the vertical.

#### (iii) Separable solution

The problem permits a self-similar solution corresponding to a separable form in polar coordinates (*r*,*ψ*) [24]. We set
*Ξ* is an undetermined length scale. The stress equations then reduce to

There are singular points in the system (3.16) and (3.17) where

For *ψ*→−*π*, we impose *Θ*(−*π*)=*θ*_{f}. For the perfectly rough case when *θ*_{f}=−*ε*, this second boundary is also a singular point. Moreover, we cannot impose regularity in view of the imposition of all the boundary conditions already. Instead, the derivative of the solution diverges and we exploit the local solution,
*π*. Numerical solutions for rough, semi-rough and smooth cylinders are shown in figure 6.

Given *ζ*→0 for *ψ*→−*π*, implying *ξ*→−1. The resulting curves (*x*,*z*)=*Ξ*(*ξ*_{α,β}(*ψ*),*ζ*_{α,β}(*ψ*)), corresponding to the two choices of sign in (3.20), generate the slip-line field through the scaling *Ξ*. Examples are shown in figure 7.

#### (iv) Collins construct

As seen below, the surface forces acting on the cylinder that are predicted by the separable solution cannot be adjusted to ensure that *F*_{x}=*O*(*ρ*_{s}*ga*^{3}/*R*) in accordance with the cylinder dynamics discussed in §3a. This motivates a modification of the separable solution that follows Collins’ solution for a cohesive ideal plastic deforming underneath a roller [26]. We attach a plug of grains to the surface of the cylinder, such that the roller and plug rigidly rotate about a point (*x*_{o},*z*_{o}) beneath the granular surface. For the back surface to depart tangentially from the cylinder, the rotation centre must lie along the line connecting the back contact point with the centre of the cylinder. If we define *Υ* to be the angle this line makes with the vertical, and *a*ℓ to be the distance from (*x*_{o},*z*_{o}) to the back contact point, then we have
*a*≪*R* and *Υ*=*O*(*a*/*R*), (*x*_{o},*z*_{o})∼−*a*(1,ℓ), and

Ahead of the plug, the material deforms plastically and the separable solution applies. Because the border between the rotating plug and this plastic region is a yield surface, it must also follow a slip line; we choose this to be a specific one of the self-similar *β*-lines intersecting the contact surface at (−*x*_{β},0) (figure 8). This *β*-line continues down to the point *P* where it intersects the lowest *α*-line, which bounds the entire yielded region. The left-hand border of the rotating plug divides that region from the immobile grains below, and for a Mohr–Coulomb material following an associated flow rule, must form a logarithmic spiral maintaining a pitch angle *ϕ* with the circular velocity of the plug material [30]. The shape of the log-spiral is given by

where *χ* is the angular position relative to (*x*_{o},*z*_{o}), measured clockwise from the vertical. This log-spiral continues the lowest *α*-line to the origin and meets that characteristic tangentially at the point *P*. This geometry demands that
*ψ*_{p} is the polar angle of the intersection point and *χ*_{p}=*ε*−*Θ*(*ψ*_{p})−*ϕ*. The lowest *α*-line is given by (*x*,*z*)=*x*_{α}(*ξ*_{α},*ζ*_{α}), and the *β*-line bordering the rotating plug is (*x*,*z*)=*x*_{β}(*ξ*_{β},*ζ*_{β}), where

The stress along the failure surface underlying the rotating plug is not given by the separable solution. Instead, the requirement that the log-spiral continues the lowest *α*-line, in conjunction with the yield criterion, implies that
*θ*=*ε*−*χ*−*ϕ*. Equation (3.26) can be integrated from *χ*=*χ*_{p} up to *χ*=0, assuming that the pressure is continuous and matches with that of the separable solution at *P* (no boundary condition is therefore applied at the back contact point, in line with previous constructions [21,22]).

#### (v) Forces and torques

We let *β*-line that border the plug with the section of the contact surface from (−*x*_{β},0) to the origin. The forces acting on this contour are the same as the forces acting on the cylinder modulo the contributions from the weight of the rotating plug. We therefore write the forces in (3.4) and (3.7) as
*A* of the plug, and where **n** denotes the outward normal to *p* and *θ*. As the pressure scales with the weight of the material, *ρ*_{s}*ga*, the contact length *a* can be scaled out and the resulting forces written in terms of dimensionless functions as
*a* for a perfectly rough cylinder with *ϕ*=36^{°}. The horizontal and vertical stresses along *x*>−*x*_{β} in (b); the redistribution of stress along the contour

As argued earlier, the angular and horizontal accelerations can only be consistent with one another provided that *F*_{x} is small. In particular, we must have _{*} which gives

To leading order we therefore find
*ϕ* in figure 9*b* for a perfectly rough and a semi-rough cylinder. Even though the semi-rough case is almost smooth (*θ*_{f}=−*ε*/10), the change in the dimensionless force function is very slight. Thus, the rolling resistance is expected to be insensitive to the roughness of the surface of the cylinder. One could also envisiage a situation in which the cylinder slides over the granular bed (and the rotating plug) in the manner of true slip when the net normal force over the contact line is small in comparison to the total tangential force. This situation is ruled out for surface friction coefficients of order unity when ℓ∼ℓ_{*}, as this guarantees that the tangential force on the surface is *O*(*a*/*R*) relative to the normal force.

In appendix A, we proceed further with the model and also compute the velocity field of the deforming granular bed, which allows one to determine the shape of the free surface in the bow wave ahead of the roller. Two significant problems emerge with this construction when using the traditional Mohr–Coulomb model: the leading-order velocity field diverges at the front contact point and the free surface descends unphysically below the bottom of the rolling cylinder. Both problems are absent for a cohesive material (appendix B), and arise here because the standard Mohr–Coulomb model predicts an excessive degree of dilation during flow, which is a known deficiency of the theory [31,32] and can be cured by introducing a relatively small angle of dilation (appendix A).

## 4. Analysis

To compare with the predictions of the plasticity theory, we fit the experimental data to model trajectories given by a force law with constant and quadratic velocity terms:
^{−1} to prevent stopping effects such as the rollback in figure 4*b* from affecting the parameters of the fit. We then performed a nonlinear least-squares fit to the solutions of (4.1) on all trials for each cylinder and granular material combination, using the initial position and velocity for each trial and the force-law coefficients *C*_{0} and *C*_{2} as parameters. Uncertainties in the fitted parameters were determined using the covariances from the nonlinear fit, scaled by the square root of the 95th percentile of the corresponding *χ*^{2} distribution, with the uncertainties in the position data all taken to be 3 mm.

The fitting results for trials on the glass sheet, sand and ballotini are shown in figure 11. The fits for the glass sheet show a small constant resistance and a more substantial quadratic component, suggesting that air drag provides the dominant resistance in these trials. In this case, *C*_{2} is related to the cylinder’s drag coefficient *C*_{d} by
*C*_{2} values suggests that *C*_{d}=2.7±0.3 for cylinders rolling on a flat surface.

For most cylinders on both the ballotini and sand, the fitted values of *C*_{0} and *C*_{2} indicate that the velocity-dependent component of the acceleration is less significant than the constant drag at the typical velocities of our trials. The ballotini shows a substantially larger constant resistance than the sand, but the two typically show comparable quadratic components. In both cases, the fitted *C*_{2} values are roughly an order of magnitude larger than those expected from air resistance, indicating that there is some velocity-dependent rolling resistance for these materials.

If the cylinder is in vertical equilibrium (*F*_{z}=*Mg*) and is fully rough, the leading-order vertical force predicted by the model can be used to solve for the contact length *a* as
*a*/*R*≈0.1 and 0.4 for our experiments with sand and ballotini, respectively. The leading-order prediction for the rolling resistance can now be written as the constant acceleration:
*Γ* in (4.4) is plotted as a function of internal friction angle in figure 12*a*. Somewhat counterintuitively, this factor decreases with the internal friction, because although the drag force *F*_{r} increases with *ϕ*, the contact area needed to support the weight of the cylinder decreases more strongly with friction. In other words, the roller digs into the medium more for lower *ϕ*, increasing the rolling resistance, as seen in the experiments where the cylinders encounter greater resistance when rolling over the ballotini in comparison to sand. A more detailed comparison of theory and experiment is shown in figure 12*b*, which shows the ratio between the fitted *C*_{0} values for sand and ballotini and the model prediction in (4.4), including the uncertainty in the friction angle for each material. The theory under-predicts the fitted resistances by roughly 50% for both the sand and the ballotini, but correctly captures the overall magnitude of the rolling resistance and its scaling with the cylinder’s mean density.

Unfortunately, the trajectory fits typically show strong covariances between *C*_{2} and the initial velocities of the faster trials, which limits our confidence in the inferred values of these parameters. Moreover, different forms for the velocity dependence of the force law (such as a linear form) lead to trajectories that do not fit the data any more or less conclusively. We cannot therefore characterize the velocity dependence of the resistance in any further detail.

## 5. Conclusion

We have experimentally investigated the resistance experienced by cylinders rolling on two dry granular media, a relatively coarse sand and ballotini (glass spheres). The cylinders were relatively long and light, with low mean densities compared to the granular materials, which distinguishes our study from some earlier experiments with spheres [17] and the bulk of the literature on wheels. Thus, surface deformations are expected to be localized and superficial. This allows us to asymptotically model the bed as being two dimensional, the deforming contact region as being much smaller than the cylinder radius, and the deflection of the free surface as being slight [19–21]. We further adopted the cohesionless Mohr–Coulomb law to describe the granular rheology. The resulting model predicts that the rolling resistance is independent of cylinder velocity, scaling with the square root of the ratio of mean densities.

The trajectories of cylinders in the experiment were fitted using a simple drag law in which the rolling resistance consisted of an overall constant plus a rate-dependent part that was quadratic in cylinder speed, similar to that expected for air drag. The measured drag was dominated by the constant term at the velocities typical of our experiments. However, the quadratic component was larger than that expected for air drag, which was estimated by rolling the cylinders down a flat glass surface. Thus, the rolling resistance appears to contain an effect stemming from rate-dependent granular rheology. We did not, however, attempt to include such effects into the rheology of our theoretical model, as might be feasible by exploiting recently proposed semi-empirical friction laws [18].

Nevertheless, the basic drag constant measured in the experiments is, to within a factor of two, consistent with that predicted by the Mohr–Coulomb plasticity model, and matches the predicted dependence on the ratio of mean densities. One important feature of the model is the characteristic flow pattern underneath the rolling cylinder. Following Collins [26], we have represented this pattern in terms of a region of plastic deformation below the front contact point which extends up to the free surface ahead of the cylinder, and back to a plug that is rigidly attached to the roller. The plug rigidly rotates with the cylinder and slides along a log-spiral failure surface over the static grains underneath. The rotating plug incorporates a finite, but small, effective slip of the cylinder as it rolls due to granular deformation, as in Reynolds’ [5] original image of the mechanism underlying rolling friction. Crucially, the plug also reorients the forces that are exerted by the bed such that the net horizontal force on the cylinder can be made small in accordance with its equations of motion. Without this reorientation, there is a physical inconsistency between the angular and horizontal momentum equations that precludes a rolling state.

Several extensions to the experiments may allow a better determination of rolling resistance of light cylinders on granular media: a longer bed would permit one to record trajectories with a larger range of velocities, and further improvement to the scene calibration would allow better inference of the cylinder trajectories in physical space, both of which would help to constrain the velocity dependence of the drag. An experimental set-up with a more precise launching mechanism and a faster way to resurface the bed would help eliminate unsuccessful trials and reduce their lateral motion and history dependence. However, such effects are interesting in their own right and may be important to any practical applications of rolling over granular media. In particular, the sideways veering of the cylinders as they leave the centre line of the channel is suggestive of an instability that deserves further exploration.

## Data accessibility

The trial metadata and object trajectories are available online at http://data.keaton-burns.com/rolling_data/.

## Authors' contributions

K.J.B. performed the experiments, performed the data processing and analysis, assisted with the modelling and helped draft the manuscript. N.J.B. supervised the experiments, assisted with the modelling and drafted the manuscript. I.J.H. assisted with the experiments, performed the modelling and helped draft the manuscript.

## Competing interests

We have no competing interests.

## Funding

Experiments and modelling were conducted at the 2016 Geophysical Fluid Dynamics Summer Study Program, Woods Hole Oceanographic Institution, which is supported by the US National Science Foundation and US Office of Naval Research.

## Acknowledgements

We thank Anders Jensen for assistance constructing the experimental apparatus. We thank Jim Hambleton for critical comments and pointing out an error in the first draft of this paper.

## Appendix A. Velocity field and free surface shape

In a generalization of Mohr–Coulomb theory, the characteristic curves of the velocity field (*u*,*w*) are given by
*ν* is the angle of dilation such that
*ν*=*ϕ*, as in the traditional Mohr–Coulomb model, do the characteristics of the stress and velocity fields coincide.

To construct the velocity fields associated with our leading-order stress field solutions, we integrate (A 1) and (A 2) numerically using centred differences on a grid defined by the characteristics of the velocity field, exploiting the self-similar solution for *θ*. We begin the integration along the lowest *α*-line within the yielded region and the *β*-line forming the yielded edge of the rotating plug. If these stress characteristics are not also velocity characteristics (*ν*≠*ϕ*), the velocity must be prescribed on both of these borders to match with the stationary unyielded material below and the rigidly rotating plug to the back. If *ν*=*ϕ*, the borders are velocity characteristics and slip is allowed along them, provided it occurs at an angle *ϕ* to the characteristics. The numerical integration then progresses upwards and rightwards to furnish the velocity field over the characteristic web underneath the last velocity characteristic leaving the top right corner of the rotating plug. Above this curve, the characteristics no longer begin from the plug, but from the yielded section of the *x*-axis. There, the *z*=0 intercepts of the upcoming characteristics must be found, and new downward characteristics launched to complete the web throughout the yielded region, taking the prescribed normal velocity *w*=−*Ω*(*a*+*x*) as the boundary condition.

The velocity field computed for the parameter values of figure 8 is shown in figure 13 for *ν*=*ϕ* and *ν*=0. In the first case, the velocity becomes arbitrarily large at the right-hand edge of the contact region, as for the problem of a critically loaded foundation [25], which reflects the well-known deficiency of the traditional Mohr–Coulomb model in predicting excessive dilation during flow. Neglecting inertial terms in the bed dynamics requires that the sand velocities are *ν*=0 furnishes an incompressible flow with velocities *ν*≠*ϕ*.

Let *z*=*h*=*O*(*a*^{2}/*R*) denote the position of the granular surface. For nearly steady rolling (*h*_{b}=*h*(−*a*), *x*_{f} denotes the front of the bow wave (where *h*=0), and we have matched to the contact arc below the roller, where *h*=*h*_{b}+(*x*+*a*)^{2}/(2*R*) for −*a*<*x*<0. Note that
*h*_{b}=0). Conversely, any dilation implies that the incoming granular surface lies below the roller (*h*_{b}>0), for which there was no experimental evidence (if anything, the medium is compacted and the cylinders roll at a position below the incoming surface). The predicted free surfaces for the two values of dilation angle are also presented in figure 13. The excessive dilation of the Mohr–Coulomb model with *ν*=*ϕ* leads to an unphysical descent of the free surface, which is avoided by taking *ν*=0.

## Appendix B. The cohesive Collins solution for *a*≪*R*

For cohesive material without gravity, the *c*, and the slip-line solution consists of a Rankine wedge leading a centred fan with circular *α*-lines and straight radial *β*-lines. The fan geometry demands that the *β*-lines be straight everywhere, and hence the rigidly rotating plug attached to the roller spans the entire contact area and *p*=*c* (*ψ*_{p} from *p*+2*cθ*=*c*(1+*π*), or *α*-line, *p*+2*cθ*=*c*(1+*π*), or _{*}≈0.43, for which the horizontal force vanishes, leading to the lift force, *F*_{z}∼4.90*ac*, and rolling resistance, *F*_{r}∼3.25*a*^{2}*c*/*R*.

Within the yielded region, the velocity field is directed along the *α*-lines with speed *du*_{α}=−*u*_{β} *dθ* and *du*_{β}=*u*_{α} *dθ* along the two characteristics, and *u*_{β}=0 along the lowest *α*-line, implying *u*_{β}=0 everywhere); cf. figure 14*b*. The velocity at the free surface is therefore *h*=0, confirming that the bottom of the roller lies at the level of the undeformed surface (figure 14).

- Received May 29, 2017.
- Accepted October 25, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.