Relaxation Times in Simple Shear and the Role of Walls

. We study the relaxation time of granular media in simple shear by means of DEM simulations (the methods being Molecular Dynamics as well as Contact Dynamics) in two and three dimensions with rough and with smooth frictional walls. While the system with rough walls behaves according to its steady state constitutive laws, the systems with smooth walls show a much stronger increase of the relaxation time with driving strength than to be expected. Employing a dynamic non-local rheology model[1] allows for a stronger increase but not su ﬃ ciently so.


Introduction
The long known non-Newtonian character of granular media [2] manifests itself in observations like a quadratic relation between stress and shear rate (Bagnold scaling [3]) or normal-stress differences [4,5]. Such quantities can readily be measured in the steady state, but of course deviations from a constant viscosity influence transients like relaxation towards the steady state, too. In this work, we're studying the relaxation times of a simple shear situation (in two and three dimensions) for different boundary conditions (rough and smooth walls) and different simulation methods: Molecular Dynamics (MD) with soft particles [6] and Contact Dynamics (CD) with rigid particles [7]. After explaining the theoretical predictions for the asymptotic relaxation time by virtue of established constitutive relations in section 1.1, we present setup and results for the two-dimensional system with rough walls in section 2 and for the systems with smooth walls in section 3. For the latter we discuss the ostensible necessity of taking into account non-local rheology in section 4 and conclude our findings in section 5.

Diffusive Evolution
Assuming a homogeneous normal stress σ zz in simple shear (driven by walls with normals in z direction), the evolution of the velocity field v(z, t) is given by with ρ being a particle's mass density. The volume fraction ν and the macroscopic friction μ = σ xz /σ zz are functions of the inertial number [8] e-mail: lothar.brendel@uni-due.de where d is the particle diameter. If ν is constant and μ(I) a linear function, (1) turns into a simple diffusion equation; otherwise it gets non-linear. But for walls with velocities ±v W having driven the system into the asymptotic regime, i.e. the "diffusion constant" d 2 /τ i just gets scaled by a function of the final shear rate. The slowest relaxation mode has wavelength 2L(I ∞ ), where the I ∞ -dependence is due to dilatancy, which also implies L(I ∞ )ν(I ∞ ) to be constant. With that we get That means, a decreasing μ and a decreasing ν both give rise to an increase of the relaxation time. For local constitutive relations [9] μ(I) = μ s + α 1/I + 1/I c and ν(I) = ν 0 − βI , (5) this factor is (1 + I ∞ /I c ) 2 /(1 − I ∞ β/ν 0 ). Considering directly the evolution of I instead of v, the counterpart to (1) is a bit more complicated, but its linearization corresponds exactly to (3), i.e. we can treat relaxation of v and of I on the same footing.

Simulation Setup
Two dimensional simulations (xz plane) in the simple shear geometry were performed using the LAMMPS[10] simulation package. As natural units we chose average particle diameter d, confining pressure σ zz and particle mass density ρ. Rough walls were made up from randomly placed disks of size 1.2 so that the effective friction coefficient of the wall was maximal [11]. The walls were kept together with constant stress in the z direction and were given shear velocities ±v W . Particles had diameters d ∈ [0.8, 1.2] with a uniform distribution, a two dimensional density of 2/3, and a microscopic friction coefficient of μ mic =0.3. The total kinetic energy in particle translation (E kin ) and rotation (E rot ) was recorded together with the shear stress σ xz . The bulk shear rate was defined to be the average over the middle third of the system.

Results
After a transient, the system reaches a steady state characterized by a linear velocity profile. The stationary shear rate is thusγ ∞ =2v W /L z , the inertial number I ∞ = √ 2/3γ ∞ . We measured constitutive relations of μ(I ∞ ) and ν(I ∞ ) which are in good agreement with the empirical predictions [9] as expressed in Eq. (5), but as shown in Fig. 1, there are alternatives like μ = μ s +αI c arctan(I ∞ /I c ), which produces a factor μ (0)/μ (  Analogous to [12] we find that all quantities reach their stationary values in an exponential way, which enables us to define characteristic times τ kin (I ∞ ) for the kinetic energy, τ rot (I ∞ ) for the rotational energy, and τ I (I ∞ ) for the bulk shear rate. Figure 2 shows the dependence of the relaxation times as function of the inertial number. The kinetic energy and the bulk inertial number behave in essentially the same way, and their relaxation time can be fairly well described by Eq. (4) and the empirical arctan-law for μ, better than according to Eq.(5). While the difference in Fig. 1 seems minute, the curvature in μ matters.
Even though the diffusion model Eq. (3) describes the system quite well, it is interesting to note that the characteristic time for the rotational energy is different from the kinetic one indicating a different underlying mechanism. In Fig. 3 E kin (t) and E rot (t) are shown for three different positions. We find that initially the rotational energy is high at the walls and low in the middle but they approach the same value. The kinetic energy is always increasing at all positions, but the stationary values are different according to the velocity profile. It seems that the rotation is responsible for the sub-linear growth of μ(I) by switching to a new shear mode with rotating particles and less friction for larger shear rates. This is further manifested in the increasing relaxation time. This rotating layer of particles close to the wall decrease the effective friction coefficient and thus the stress in the bulk making the diffusion process slower as shown in Fig. 2.

Simulation Setup
Smooth walls represent a stronger perturbation of local packing geometry, and they can transmit force and torque only via contact friction: Sliding friction is set to μ mic = 0.5 everywhere, rolling and torsion friction are absent.
In the three-dimensional Contact Dynamics simulations, N = 15300 particles with diameters ∈ [0.8, 1.0] and zero coefficient of restitution are placed at rest in a cuboidal box with a base of size L x × L y = 10 × 10 and periodic boundary conditions. In z-direction, walls apply the unit pressure σ zz by means of a constant force and have a fluctuating distance of L z ≈ 100 and a fixed velocity ±v W in x-direction; cf. [12] for further details. We also carried out 2D MD simulations with smooth walls, the parameters of which correspond essentially to the ones in the previous section, except for μ mic = 0.5 at contacts between particles and to the walls.

Results
In contrast to walls made out of particles, the smooth walls allow for a finite slip velocity, which itself is a dynamical quantity and is a non-trivial response to the driving velocity. Nevertheless, well defined boundary conditions for the bulk evolution are established in the sense that about 9 particle diameters away from each wall (defining an effective system size L z < L z , which we will denote L in the following), the shear rateγ and thus the inertial number takes on a constant value I 0 (apart from fluctuations, of course) after less than 100 natural time units (not shown). The same time scale actually applies to the total rotational energy. During the subsequent relaxation of the bulk towards this shear rate, the corresponding velocity difference  grows just as the slip velocity decreases. Due to the wall friction being Coulombian, this decrease has no influence on the driving force . Brief, the finite slip velocity provides a "buffer" as to enable I 0 as boundary condition (and thus Neumann BCs for the velocity).
In the steady state, the constitutive relations are measured, as shown in Fig. 4. Again, we measured the relaxation time of the bulk inertial number, the results of which is displayed in Fig. 5, for the 3D CD simulations as well as for the 2D MD simulations.
As already mentioned in section 1.1, a linear μ(I) does not contribute to an I ∞ -dependency of the relaxation time, leaving only the very weak influence of ν(I), which, as shown in Fig. 5, is essentially negligible. A logarithmic law μ(I) = μ s + ln(1 + αI), on the other hand, actually leads to a factor μ (0)/μ (I ∞ ) = 1 + αI ∞ , which has the correct functional form seen in the simulation data in Fig. 5, but the slope is too small by a factor of six.
We omitted the comparison for the two-dimensional MD-simulations, because the discrepancy is the same. But we include the data to demonstrate that the strong increase of τ with I ∞ is quite generic and no artifact due to rigid particles or perfect collisional dissipation.  The 3D data is compared to the prediction according to Eq. (4) using the fits, linear and log, from Fig. 4. The green curve corresponds to a linear μ(I ∞ ) with an increase of τ solely due to dilatancy. The blue curve incorporates the effect of a negative curvature in μ(I ∞ ), but the slope cannot match the simulation data. The normalization factor τ 0 is τ(0) from Eq. (4).
Even though we consider boundary conditions for I given within the system, about 9 particle diameters away from the actual walls, we cannot exclude any additional influence of the walls. They cause inhomogeneities, which have been shown to be captured successfully by non-local rheology [13].

Non-local Rheology
In the time dependent non-local rheology according to [1], an additional field g, the granular fluidity, is introduced to describe the shear response I = τ iγ = g μ 1 and obeys where the parameters t 0 and A set the intrinsic time and length scale of the evolution of g. For steady, homogeneous flow (constant g > 0) it yields g = (μ − μ s )/(αμ) which corresponds to I c → ∞ in (5). We can eliminate the field μ from the equations to keep only I and g: Then, linearizing around I ∞ and g ∞ = I ∞ /(μ s + αI ∞ ) ≡ I ∞ /μ ∞ leaves us with where ν ∞ = ν(I ∞ ). A solution with vanishing δI and δg at the boundaries (i.e. Dirichlet BCs for I and g) is and yields (besides the irrelevant μ * ) The static limit τ nl = t 0 L 2 /(π 2 A 2 ) is peculiar, because it does not contain parameters of the constitutive relations μ(I) and ν(I), but it is non-zero as opposed to the case considered in [12], using the model from [14]. On the other hand, τ nl is bounded by τ , which is clearly exceeded by the simulation data, as shown in Fig. 6. Hence, taking into account non-local rheology allows for a larger variability of the relaxation time, but is not enough.  Figure 6. The relaxation for the 3D system with smooth walls and as predicted by the non-local rheology for L 0 = 81.5, A = 4 and t 0 = 0.6A 2 /(ν 0 α), exhibiting a non-zero τ nl (I ∞ → 0) and a pronounced, A-dependent slope for small I ∞ . But it is bounded by the result for a locally valid, linear μ(I) (cf. Fig. 5), which prevents it from agreeing with the simulation data.

Conclusion
We have studied the very simple concept of the time scale needed by a granular system in simple shear to relax into its steady state of homogeneous bulk shear rate (and thus inertial number I ∞ ). The theory, linearized momentum conservation closed by constitutive relations for macroscopic friction and volume fraction, predicts the somewhat counter-intuitive increase of the relaxation time τ with increasing driving velocity. These predictions are pretty accurate for rough walls, but are much smaller than what we measured in the simulations with smooth walls.
Employing a non-local rheology model [1] provided a stronger variation of the relaxation time compared to the local rheology, but only in one direction (lower times). The role of the boundary conditions for the granular fluidity remains to be clarified, though (Dirichlet BCs have been chosen).