Relaxation processes after instantaneous shear-rate reversal in a dense granular flow

A numerical model of granular material at different packing fractions and under steady shear is submitted to a sudden shear reversal. We monitor consequences of the strong density and shear rate spatio-temporal heterogeneities, on the constitutive relations. We show that the dynamics can be decomposed into two subsequent regimes spanning a time scale inversely proportional to the shear rate. In the first regime, the non-local constitutive relation proposed by Bouzid et al. (Phys. Rev. Lett., 111 (2013) 238301) is satisfied hence accounting for the spatial heterogeneity of the fluidity parameter. However at later time, we find that the local constitutive relation can be applied, in spite of the fact that the fluidity parameter remains heterogeneous. This result opens a discussion on the significance of fluidity and hidden variables accounting for non-locality in a frictional packing of grains.

Introduction. -The understanding of dense granular flows received an important impulse by the discovery of the local rheology laws that systematized the results of many experiments [1][2][3]. Dry granular media in avalanches or planar Couette geometries are standard benchmarks where the normal pressure P and shear stress τ can be measured and analysed separately. It was found that in hard grain systems, the normalized shear stress μ = τ P -which can be associated to an effective friction coefficient-depends only on the dimensionless inertial number , for a flow of particles of diameter d and mass density ρ p , shear rateγ and confining pressure P . The inertial number results to be the ratio between the microscopic time t 1 = d √ P/ρp to displace one grain in one diameter under pressure and the time t 2 = 1 γ associated with the macroscopic flow scale. The relation μ(I) constitutes a local rheology relation [1][2][3]. Experiments and simulations show that μ is monotonously increasing with I, with a finite limiting value for vanishing flow rate lim I→±0 μ(I) = ±μ 0 . At steady state, it is associated with a second relation: φ(I) expressing the decrease of packing fraction with the inertial number. The limit I → ±0 yields the maximum packing fraction value φ * depending on the microscopic granular details, particularly on the friction coefficients. For spherical grains, φ * would essentially decrease when granular friction increases. Note that the φ * value is obtained by extrapolation of numerical results and always implies steady and homogeneous fields as far as shear rates and volume fractions are concerned. It is important to keep in mind that the only true limit, where a solid phase would emerge as a real rigidity transition, is obtained for zero friction and is called the jamming point, sometimes called the random close packing limit [4]. Therefore, as we will see later, nothing prevents the observation, as a transient, of a granular packing flowing at packing fractions larger than φ * as long as packing fraction remains smaller than the random close packing limit.
The local rheology picture was extended to non-planar geometries by assuming that stress and shear-rate tensors remain parallel with a proportionality constant associated to μ, with I computed from the invariants of the tensor [5]. Computer simulations of a rotating drum show that those tensors are indeed not parallel but still their invariants are related via the μ(I) rheology [6]. Simulations of a flow on an inclined plane have shown also a small but finite angular deviation when an objective analysis of the tensors is performed [7]. However, for non-homogeneous fields, non-local effects were observed 64002-p1 either close to jamming [1,8,9] or in the vicinity of a shear band [10][11][12]. Different theoretical visions were proposed to model or interpret these features [9,[13][14][15][16]. In the approach by Bouzid et al. [16], the local rheology relations were modified to take into account these effects in the form of a non-local rheological law μ = μ(I) [1 − χ(κ)], where κ = d 2 ∇ 2 I/I. Using this model, it has been shown that a relaxation length develops from the shear band, which diverges at the maximum packing fraction φ * [16,17].
In the application of the rheological laws (local and nonlocal) to model flow regimes and hydrodynamic fields it is assumed that they would retain their form under timedependent conditions and, therefore, could be used as input in the momentum conservation equations. Here, we investigate if this hypothesis is valid in a transient regime.
With this objective in mind we analyse the rheological response in time-dependent planar Couette flows using molecular dynamics simulations of a frictional granular packing. It is found that the system evolves in two distinct phases. In the first one the non-local rheology accurately describes the system evolution, with similar values for the non-local correction χ(κ) as compared to refs. [16,17]. In the second phase the temporal evolution is characterised by the local rheology, even though κ does not vanish.
Simulation method. -The numerical set-up consists in a two-dimensional dry granular medium confined between two rough walls in the absence of gravity (see fig. 1(a)). The medium is made of a polydisperse mixture of circular grains with uniform mass density and diameters that are uniformly distributed in the range [0.5d, 1.5d], where d is the average diameter and m p is the average grain mass. Such a mixture allows us to have a stable amorphous structure that does not present segregation nor crystallization. The walls are made of similar particles, which will be forced to move at imposed velocities ±U w to produce a planar Couette flow. The shear rate iṡ γ w = 2Uw H , H being the height of the flowing region. Transient regimes are obtained by instantly reversing the imposed velocities on the walls (see fig. 1(b)). The rapid reversion of the shear rate generates a strong perturbation that generates important inhomogeneities in the relevant fields when relaxing to their stationary value. We opted for this perturbation instead of a softer one (for example, using a sinusoidal oscillation) because it does not introduce additional time scales. Here, we will analyse the relaxation process in the dimensionless timeγ w t, with t measured from the shear-rate reversion.
The model for the grains is the one used in ref. [18], which we summarize here briefly. The grains interact with visco-elastic forces in the normal and tangential directions, added to a tangential Coulomb friction. The normal and tangential elastic forces are linear with elastic constants k n and k t , respectively. In the limit of high elastic constants, the microscopic dynamics is characterized by a unique time scale given by the shear rate. Therefore, the intensity of the shear rate can be fixed to set units. The control parameter is the volume fraction φ which is set by a protocol described below. The spring constant for the normal force is set to k n = 10 10 m pγ 2 w , which is much larger than the maximal pressure reached in the simulations and therefore we remain in the hard particle limit for which the μ(I) rheology was obtained. The tangential spring constant is fixed to k t = 0.5k n and the dissipative terms are fixed to obtain restitution coefficient in the range [0.1, 0.9]. Finally, the Coulomb friction coefficient is chosen equal to μ p = 0.4. The equations are integrated using a time step that is 1/25-th of the average collision time, small enough to guarantee that particle collisions are correctly sampled.
The system is periodic in the horizontal x-direction, while it is limited by the walls in the vertical one. The horizontal size is L = 45d and we studied two system heights: H 1 = 18d and H 2 = 40d. Finally, the wall thickness is W = 1.5d, such that all particles with their center in this region move together with the wall.
The initial configurations are prepared as follows. N particles are placed in the system between the walls (N = 800 for H 1 and N = 1800 for H 2 ). A fixed vertical stress P e is imposed on the walls to compress the system and a shearing motion is added until the system relaxes to the steady state. Depending on the imposed pressure, different global volume fractions φ G between the walls are obtained in the range 0.767 < φ G < 0.805 for H 1 and 0.778 < φ G < 0.809 for H 2 .
Having obtained the initial configurations, the wall separation is fixed and the following simulations are done at constant volume. Planar Couette flows are simulated with an imposed wall velocity U w , until a stationary state is achieved. The stationary state is slightly inhomogeneous in the vertical direction (see fig. 2). Hence, we measure in the central region of the box (averaged over one third of the system height) the steady-state volume fraction φ c , friction coefficient μ c , shear rateγ c , and inertial number I c . The friction coefficient depends, as usual, on the inertial number and the results are fitted to the standard expression μ c = μ 1 + μ2−μ1 1+I0/Ic [1], where the fit parameters are μ 1 = 0.277, μ 2 = 0.85 and I 0 = 0.364, the same as those obtained in ref. [19]. Also, the measured vol- For each cell width it is then possible to extrapolate the global maximum packing fraction values φ * G = 0.805 for H 1 and φ * G = 0.813 for H 2 . We proceed then to produce a planar Couette flow with wall velocity ±U w , where the sign alternates periodically ( fig. 1(b)). The alternation period is large allowing the system to reach a steady state before it is perturbed by the sudden change in shear direction. The process is repeated and averaged considering 1000 cycles for H 1 and 250 cycles for H 2 . The purpose of this study is to characterize from the rheological point of view, the relaxation process to the new state after a sudden change of the shear-rate sign.
System evolution. -The relaxation process is strongly inhomogeneous temporally and spatially, however it can be decomposed into two distinct phases. Figure 2 displays the temporal evolution of the volume fraction φ(y, t) (figs. 2(a) and (b)) and the velocity u(y, t) profiles ( fig. 2(c)). All profiles are obtained as averages over spatial bins in y and over the cycles of shear-rate reversions. The inertial number profiles are obtained by locally fitting the velocity profile to low-degree polynomials, to later derive them.
Throughout all the relaxation numerical experiments we performed, we noticed that the system evolution can be characterized by two subsequent time scales which can To highlight the evolution in the two phases, for illustration purposes only, in the top figure we offset on both axes the relaxation curve for φG = 0.776. First, during τ1 the system does not follow the local rheology and the volume fraction increases until a maximum value is reached. In the second phase, of duration τ2, the local rheology μ(I) is fulfilled but the volume fraction, which decreases monotonically while I increases, is not given by the steady-state relation φ(I).
be simply visualized by inspection of the density relaxation profiles. First, just after the shear reversal and due to the increase of pressure at the wall, a transverse flux of 64002-p3 matter increases the density towards the centre. Then the density in the centre reaches a maximum and thereafter, in the second phase, start to relax down to its steadystate value. Both phases scale with 1/γ, however, the second phase lasts about ten times longer than the first phase. After the first phase, the pressure field becomes spatially homogeneous and reaches its steady-state value (not shown). Thus, the density relaxation is not driven anymore by pressure gradients but qualitatively, by processes more akin to a diffusional relaxation, associated to non-affine motion. At all times, the corresponding shearrate fieldsγ(y, t) are spatially heterogeneous (see fig. 2(c)).
Local rheology. -Now we investigate the relations linking the dynamical friction coefficient and the packing fraction, to the inertial number. For this purpose we measure the volume fraction φ c (t), inertial number I c (t), and friction coefficient μ c (t) in the central region of the simulation box (averaged over one third of the system height, namely in −3d ≤ y ≤ 3d for H 1 and in −6.7d ≤ y ≤ 6.7d for H 2 ). Figure 3 presents the system evolution of those variables in the μ-I and the φ-I spaces, for different steady-state volume fractions φ G . As a reference, in black, we display the relations μ(I) ( fig. 3(a)) and φ(I) ( fig. 3(b)) obtained in the central region for the stationary state. After the velocity reversal, the first phase corresponds to a strong departure from the μ(I) relation. However in the second phase, the relaxation to steady state is taking place on the μ(I) curve. This final relaxation can be done either by increase or decrease of I c according to the value of the global packing fraction φ G .
For the packing fraction evolution (see fig. 3(b)), in no instance of the relaxation process, the local φ(I) relation is satisfied. Interestingly, the maximal φ c value reached in the centre is always higher than the steady-state value φ * but still, I c remains finite. However, in all cases, the volume fractions are smaller than the random close packing limit.
Non-local rheology. -To measure quantitatively the deviation from the stationary rheology, we define in the centre the parameter which is presented in fig. 4(a) as a function of time for φ G = 0.767. To test the non-local rheology in the spirit of Bouzid et al. [16], we computed the parameter κ = d 2 ∇ 2 I/I associated with the local fluidity relative to the neighborhood. We obtained in practice a good value of the second derivative of I in the centre, by fitting I(y) in a central region of height 12d by an even 4th-degree polynomial. Figure 4(b) displays for φ G = 0.767 the time evolution of this parameter and one sees that it reaches zero after a time τ 1 but thereafter displays a rebound before reaching zero in the second part of the relaxation dynamics. Remarkably, when Δμ is displayed as a function of κ = d 2 ∇ 2 I/I (see fig. 4(c)), one obtains a good collapse for the different volume fractions on a single curve χ(κ), hence validating the non-local relation as proposed by Bouzid et al. [16] of the type Note that here the relation goes beyond the linear limit (Δμ ∝ κ) as presented in ref. [16] but it compares well with the extended data obtained for a sheared stationary flows in [17]. However for t > τ 1 , the non-local extension to the rheology is not valid. The departure in κ from zero is small but, nevertheless, systematic.
Time scales. -We have shown that the system relaxation takes place in two phases. Note that the corresponding time scales reflect complicated processes involving the whole systems dynamics. In this report we do not seek to solve the complete dynamics but only to characterize the empirical behaviour and in particular the phenomenology around the point where the system reaches maximum packing fraction. In the first phase, the density increases in the centre and the departure to locality Δμ(t) goes to zero. When plotted Δμ(t) in log-linear scale (not shown), it decays exponentially for all volume fractions, allowing us to measure the duration of the first phase τ 1 by an exponential fit.
Note that the time to reach the maximal value is related to τ 1 and for φ G close to φ * G , ones obtained a time scale about 3 times larger (not shown). The second time scale τ 2 is obtained by systematically fitting the exponential tail relaxation of φ c (t) − φ c after the maximum. The initial time scale τ 1 is presented in fig. 5(a) as a function of global volume fraction φ * G − φ G and τ 2 is presented in fig. 5(b) as a function of (φ * G − φ G ) 3/2 . These time scales generally depend on the finite height H, however importantly, close to the maximum density φ * G , they remain finite and thus we do not find the signature of a critical phenomenon corresponding to a critical slowing down of the dynamics.
Summary and discussion. -The relaxation process of a granular system in a planar Couette geometry at fixed volume is studied after the imposed shear rate is instantly reversed. We found that the system evolves in two phases. First, the steady local relation μ(I) is not satisfied but the measured value approaches exponentially the steady relation, with a relaxation time that depends on the global density. In this phase the non-local rheology is satisfied with similar values for the non-local correction χ(κ) as compared to refs. [16,17]. Because of the pressure increase at the boundary there is a mass flux towards the centre of the cell such that density increases but the local constitutive relation φ(I) obtained for a steady uniform regime is not satisfied. To better understand this phase one would need to consider the motion in the vertical direction. Thereafter, there is a second phase where the system evolves following the local relation μ(I) found for steady flows. Finally, the time scales for both phases depend on the distance to the maximum packing fraction, but do not show any critical behaviour.
The salient result of this study is that the fluidity parameter κ, associated to the Laplacian of the inertial number, is a good predictor for the non-local rheology observed in the first phase, where these effects are quite important. However, in the second phase, where we identified smaller but systematic non-zero values for κ, we observe that the non-locality does not manifest itself any more in the μ(I) rheology. This observation points on the possibility of a relevant hidden variable accounting for non-locality. This variable would be similar to the fluidity parameter as defined by Bouzid et al. [16], but would relax faster to reach zero in the second phase. In some sense, the recent work of Bouzid et al. [17] on a micro-rheological test for the nonlocal rheology at steady state, reaches similar conclusions. Both works stress on the need to provide a deeper interpretation for the micro-mechanical processes at play yielding non-local behaviour for dense granular packing under shear. More research is needed to consistently include the two phases found in this work in a unified rheology that can be used as input in the momentum conservation equations for dense granular flows. * * * This research is supported by Fondecyt Grants No. 1120775 and 1140778 and Ecos-Conicyt Grant C11E04. ER acknowledges the support of a Becas CONI-CYT. BA is supported by Institut Universitaire de France. This work is funded by the ANR JamVibe and a CNES research grant.