Wall roughness and nonlinear velocity profiles in granular shear flows

Inhomogeneous velocity profiles in granular flows are well known from both experiments and simulations, and considered as a hallmark of nonlocal behavior. By means of extensive contact dynamics simulations, we show that the sigmoidal velocity profiles in 2D flows of rigid disks are controlled by the roughness of driving boundary walls. We find that the velocity profile becomes linear for a critical value of wall roughness up to an exponential decay close to the walls with a characteristic length that does not depend on the flow thickness and rate. We describe the velocity profiles by introducing a state parameter that carries wall perturbation. By assuming that the local shear rate is a linear function of the state parameter, we obtain an analytical expression that fits velocity profiles. In this model, the nonlinear velocity profiles are explained in terms of the effects of wall roughness as boundary condition for the state parameter.


Introduction
Many natural processes, such as debris flows [1,2], and industrial operations in powder technology [3], involve granular flows with various particle properties, boundary geometries and driving mechanisms.Dissipative contact interactions and the absence of long-range cohesive forces lead to a complex dynamics, in which the energy input and stresses from the boundary are balanced by internal dissipation and momentum exchange between particles.In the classical Mohr-Coulomb model, only quasi-static deformations and limit states are considered and the behavior is modeled in terms of a yield surface characterized by an isotropic effective (or internal) friction coefficient μ = σ t /σ n , where σ t is the shear stress and σ n is the normal stress [4].This model was later supplemented in soil mechanics with a flow rule by considering dilatancy as a function of a state variable measuring the distance to the critical state (steady isochor flow) and a phenomenological evolution rule of the state variable [5].
More recently, the frictional behavior was extended to inertial states, in which both the effective friction coefficient μ and packing friction Φ are described as functions of the inertial number I = γd s (ρ s /σ n ) 1/2 , where γ is the shear rate, d s is the mean particle diameter and ρ s is particle density [6,7].
A key pending issue is the observation of flow inhomogeneities, which can not be predicted or analyzed within the afore-mentioned models of granular flow.The most common example is spontaneous strain localization in shear bands of small thickness [8].In the Mohr-Coulomb framework, this phenomenon is described as a split of the e-mail: paul.schuhmacher@umontpellier.fr material into two rigid blocks slipping past each other and explained as material instability [9].In simple shear flows, it is observed that, depending on the wall roughness and shear rate, the flow may consist of a solidlike region that does not flow or creeps and a fluidlike region that flows, with an interface that may evolve in time [10].
Even when the flow occurs in the whole bulk of the material, the shear rate is often observed to be nonuniform despite uniform stress across the flow [11].This feature is obviously in contrast with the μ(I) rheology, in which μ is assumed to be a unique function of I so that a constant value of μ across the flow implies a uniform value of I and thus a uniform shear rate.Similar inhomogeneous velocity profiles have been observed in other soft glassy materials such as emulsions and colloidal suspensions [12,13].Recently, several models based on a nonlocal formulation of the rheology have been introduced to account for these observations [14][15][16].
In this paper, we use contact dynamics simulations to investigate the velocity profiles in a systematic way in stress-homogeneous 2D shear flows driven by moving rough walls and confined by a constant pressure perpendicular to the flow.We are interested in the effects of wall roughness and flow thickness.As we shall see, independently of the inertial number, wall roughness controls the velocity profile as long as the flow thickness is below 130d, where d is mean particle size.Above this thickness, the shear strain is uniform across the flow.We show that this observation can be rationalized in a model accounting for a state parameter as carrier of perturbation introduced by wall roughness.

Numerical procedures
An assembly of 6000 disks is confined between two rough horizontal walls as shown in Fig. 1.The particle size distribution P(d) is uniform in particle volumes (P(d) ∼ d −2 ) with diameters in the range [d m , d M ] where d M = 2d m .This particle size polydispersity is necessary in 2D to favor long-range geometrical disorder.The system is periodic along the flow.The walls are made rough by attaching a periodic array of particles ('wall particles') of diameter d w which can not rotate.The wall roughness is controlled by the size ratio R = d w /d s and the distance w between wall particles; see Fig. 2. The top and bottom walls are moved horizontally at constant velocities v and −v, respectively, in the horizontal direction x.The bottom wall is immobile along the vertical direction y whereas the top wall is free to move vertically and subjected to a constant compressive stress σ n = σ yy .The nominal inertial number I e defined from boundary conditions is set to 0.01.In these conditions, we are sure to avoid the intermittent regime at low inertial number in which averages and velocity profiles depend on the integration time [10].
We assume that the particles are very stiff with a Young modulus well above the confining stress σ n .In this limit, the particles may be considered as fully undeformable and simulated by the contact dynamics method (CDM), which is a discrete element method based on nonsmooth contact laws [17,18].In this framework, the equations of dynamics for all degrees of freedom are integrated together with the kinematic constraints arising from unilateral contacts and Coulomb friction law.Hence, unlike molecular dynamics method, which employs particle overlaps as strain variable, in CDM the forces are calculated simultaneously with particle velocities via an iterative implicit scheme.In our simulations, we set the friction coefficient between particles to μ s = 0.3 and the normal and tangential restitution coefficients to zero.In dense granular flows, the flow behavior is insensitive to the values of restitution coefficients since energy is mostly dissipated by frictional slipping events, and the collisions have a multi-contact nature.

Velocity profiles
In order to extract average velocity profiles from the simulation data, we map the information carried by the particles to an Eulerian regular grid of 128 × 64 points.At each point r of the grid at a given time t, the local mass and momentum are respectively given by ρ m (r, t) = i m i δ(r−r i ) and p m (r, t) = i m i v i δ(r − r i ), where r i (t) is the position of particle i at time t, m i is its mass and v i is its velocity.Equivalent continuum fields are defined by convolution with a Gaussian distribution: is a Gaussian function of width w = 2d, defining the spatial resolution.The velocity field is then given by v(r, t) = p m / ρ, which satisfies the mass and momentum conservation [19].This convolution leads to smooth derivable fields.
At a given value of the nominal inertial number I e defined from the nominal shear rate γe = 2v/h, the average height h of the flow remains constant.We divided the flow into 50 horizontal layers of thickness δh d.The velocity profiles were obtained by averaging over time along each layer.The profiles shown below were calculated in the steady state over a cumulative shear strain higher than 6.
Figure 3 shows velocity profiles for I e 0.01 and several values of the roughness angle δ w defined by (see Fig. 2): The velocity profiles clearly depend on the roughness angle and they present a sigmoidal shape, in agreement with measurements reported by several authors [11,20] .But we also see that for a critical value of this angle δ c w 0.75 radians or nearly 45 • , the profile is remarkably linear.We also have nearly the same profile for two very different values of δ w above and below δ c w .The use of roughness angle is motivated by the observation that the velocity profiles for different combinations of size ratio R and w are controlled only by the roughness angle, which combines these parameters (not shown here).It is a physically plausible parameter as the effect of roughness depends on how momentum is transmitted from wall particles to the mobile particles of the sheared sample.The amplitude is clearly given by m s v, but its effect depends on its direction.For given d w , the lowest value of δ w is d w /(d w + d s ), which vanishes only when d w /d s → 0. Figure 4 displays the shear rate γ * at the center of the flow normalized by the nominal shear rate γe for a set of simulations in which both w and d w have been changed.At the minimum value of δ w , the velocity profile has its largest deviation from the linear profile, corresponding to the lowest value of γ * .As δ w increases from its minimum value, γ * increases and reaches its highest value at δ c w .At this point, the velocity profile is linear and the shear rate at the center of flow coincides with the nominal shear rate γe .When δ w is further increased, γ * again declines.This non monotonous dependence of the velocity profile on roughness angle reveals the nature of the wall-flow interaction.Indeed, the angle δ c w 45 • corresponds to the major principal direction of the strain-rate tensor or the direction of maximal contraction.Hence, at δ w = δ c w , the direction along which the wall momenta are on the average transmitted to the granular material is consistent with that of the flow.For larger and smaller values of δ w , the mismatch between the two directions leads to lower transmission of wall momenta and hence lower shear rate at the center of the flow.
Note that in many experiments and numerical simulations, the walls are made rough by attaching contiguous particles of the same size as flowing particles to the wall.For such a 'natural' roughness, we have R = 1 and w = d w , so that the roughness angle is 30 • , which is below δ c w .According to our data, this angle leads to a nonlinear profile with γ * /γ e 0.6, thus the counterintuitive conclusion that 'natural' roughness is not the best way of setting a granular material in motion.Given the value of the critical roughness angle, for contiguous wall particles the best size ratio is R 2.3 in 2D, for which a linear velocity profile is expected.

Perturbation model
The picture briefly presented above is all the more plausible that the nonlinear velocity profiles find here their origin in the diffusion of momenta from the boundary layer.Put differently, the nonlinear profiles are a consequence of the perturbation of a linear profile by the wall and its roughness.Hence, for large enough values of flow thickness h, the velocity profile is expected to become linear.Fig. 5 displays the velocity profiles at constant values of wall roughness and inertial number but for three different values of h.We see that the profile is linear for h/d = 131.The distance over which the perturbation produced by wall roughness propagates (the exponential falloff of velocity near the walls) is essentially the same in these simulations (note that the distances are normalized by h).Most experimental or numerical studies in plane-strain geometry have been performed for a flow thickness below 40d.Thick flows (100d) were simulated by da Cruz et al. but without focusing on boundary effects [7].
In order to capture the wall effect on velocity profiles, let us consider a scalar state parameter χ that carries the effect of wall roughness across the flow.This parameter can be a descriptor of the microstructure or dynamic heterogeneities of the flow.The flow behavior will depend on both the shear rate and this state parameter.This parameter reflects the effects of the two walls so that χ(y) at a position y is the sum two contributions f (h/2 + y) and f (h/2 − y) by the same function f but depending on the distances h/2 + y and h/2 − y from the two walls: It is easy to see that all even derivatives of χ with respect to y have the same structure as in ( 2) so that all even derivatives are proportional to χ(y).The only regular function satisfying this condition is hyperbolic cosine: where χ * = χ(0) and ξ is a characteristic length, which can be interpreted as the distance over which wall perturbation propagates.It should be intrinsic, proportional to d and independent of h.The shear strain field γ(y) is assumed to depend on the state parameter.As γ(y) has the same symmetry as χ(y), it can be expanded to first order as γ(y) B + Aχ(y).The constant B can be determined as a function of shear rate γ * at the center of flow: By integrating with respect to y and given the symmetry of the velocity field, we get This functional form fits excellently our data as shown in Fig. 3.The effect of flow thickness h can be analyzed by considering equation ( 5) at y = h/2 where v This relation suggests that, for a given driving velocity v, A declines as h increases.Hence wall perturbation declines everywhere in the flow; see Fig. 5.This result is independent of the precise mechanism of wall perturbation or the state variable that carries the perturbation.It only requires that ξ be independent of system size.By analyzing the simulation data for different values of the inertial number (up to I e = 0.08), wall roughness δ w and flow thickness h, we find that ξ declines only slightly from 7d to 6d as I e increases.For a given flow thickness, on the other hand, the amplitude A declines as δ w increases and nearly vanishes at δ c w , then increases again for larger values of δ w .This is consistent with linear profiles observed at critical roughness.

Conclusion
In this paper, we used extensive contact dynamics simulations to analyze the effect of wall roughness on simple shear flows of granular materials in the steady state.It was shown that the velocity profile is linear for a flow thickness of the order of 130d.For smaller thickness, the velocity profile has a sigmoidal shape that varies with wall roughness.'Roughness angle' that encodes the geometry of wall roughness and sensitively describes the velocity profile.We showed that, independently of flow thickness, the velocity profile becomes linear at a critical value of roughness angle.Physically, we argued that this angle coincides with the most contractive direction in the flow, suggesting thus that for this critical angle momentum is optimally transmitted from the walls to mobile particles of the flow.We also introduced a simple model based on a state parameter that carries wall perturbation across the flow.Assuming linear perturbation, we found an expression for velocity profiles that accurately fits the data and explains some of the observed features.This framework suggests that a state variable can be used to enrich a purely local rheological law.
We are grateful to Jean-Noel Roux for useful discussions and references that he brought to our attention.

Figure 1 .
Figure 1.Boundary conditions for simple shear flow.The arrows represent the particle velocities.The boundary walls and particle velocities are periodic in the horizontal direction.

Figure 2 .
Figure 2. Illustration of wall roughness.The gray particles are attached to the wall while the white particle is mobile and belongs to the flow.

Figure 3 .Figure 4 .
Figure 3. Velocity profiles for four different values of roughness angle δ w and a constant nominal inertial number I e 0.01 in the steady state, for a fixed average flow thickness h/d ≈ 50.Solid lines are fits by equation (6).

Figure 5 .
Figure 5. Velocity profiles for three different values of flow thickness in the steady state at constant nominal shear rate γe and wall roughness δ w = 0.44. )