Transients in pressure-imposed shearing of dense granular suspensions

Granular materials whether dry or immersed in fluid show dilation or compaction depending upon the initial conditions, solid fraction and normal stress. Here we probe the transient response of a dense granular suspension subjected to change of applied normal stress under simple shear. In this aim, normal-stress-imposed discrete element particle simulations are developed considering the contributions arising from the drag induced on the particles by fluid phase. These pressure-imposed simulations show transient behaviors of dense granular suspensions such as dilation or compaction before reaching a steady state following the μ(J) rheology. Less expectedly, the transient behavior, in particular the height of the system as a function of applied strain, can also be described by assuming that the system follows the steady μ(J) rheology at all times.


Introduction
Understanding and predicting the rheology of dense granular suspensions is highly important to understand the behavior of e.g. debris or mud flows [1]. Additionally, it also finds application in civil, mining and food engineering sectors while handling concrete, drilling muds, petroleum extraction and manufacturing of chocolates and other food products [2]. In the last two decades, there has been substantial advancement in understanding and modeling the steady flow of granular suspensions which has led to the development of µ(J) rheology [3,4]. For an homogeneous mixture of non-Brownian and neutrally-buoyant rigid particles in a viscous fluid, dimensional analysis shows that in simple shear the rheology is described by a couple of relations expressing the macroscopic friction µ = τ/P p (with τ and P p respectively the shear stress and the particle pressure in the gradient direction) and the solid fraction φ as a function of the so-called viscous number J = η fγ /P p .
Because φ SS (J) < 0, the µ(J) rheology captures the change of volume for the particle phase between two steady states, that is, dilation (resp. compaction) under increase (resp. decrease) of J [4], a phenomenon observed in immersed soils [5,6] or submarine avalanches [1,7]. However, the transient behavior of the system during dilation or compaction is less characterized, and it is not known if an extension of the µ(J) rheology to transient states is needed to capture it. When the system dilates, the fluid is sucked into the granular packing as a result of newly created pore spaces. In contrast, during compaction, * e-mail: shivakumar.athani@univ-grenoble-alpes.fr * * e-mail: romain.mari@univ-grenoble-alpes.fr the fluid is expelled out of the void spaces and accordingly the solid fraction increases [2,8]. The associated mechanical coupling between the fluid and particle phases happen before reaching the steady state where the µ(J) rheology applies.
In this work we start filling this gap with numerical simulations of a sheared suspension undergoing dilation or compaction. While there are many numerical studies regarding the transient response of a sheared granular suspension, all of them address the problem under fixed volume conditions [9][10][11][12]. Conversely, though several experiments [4,[13][14][15] and particle simulations based on discrete element method (DEM) [16][17][18][19][20] have employed pressure-imposed conditions, they have been restricted to steady state investigations. We therefore developed DEM simulations of a simple shear of a granular suspension under imposed external particle normal stress P ext . One way to capture the transient behavior is to impose a step-change in J ext ≡ η fγ P ext , which leads to either dilation or compaction of the suspension. We show that in this setup, the entire dynamics of volumetric change can be well described by the µ(J) rheology, even though this rheology is a priori only valid in steady state.

Numerical method
Using a Discrete Element Method (DEM), we simulate a monolayer system of N = 1000 spherical particles immersed in a Newtonian fluid with viscosity η f , subjected to an homogeneous simple shear at rateγ between two horizontal walls made of frozen particles, as depicted in Fig. 1. The choice of the monolayer, instead of a threedimensional system, gives us access to larger separations between the two walls, which limits finite-size effects. Dimensions of the system are L x and h in the horizontal and vertical directions, respectively, and the simulation Figure 1. Our 2D monolayer setup. The system (grey particles) is sheared at fixed applied vertical pressure P ext between an impermeable fixed bottom wall (blue particles) and a permeable top wall free to move vertically to adjust the system area.
domain is periodic in the horizontal direction. Particles are bidisperse, with small particles of radius a and large ones of radius 1.4a, mixed in equal area proportion. Because of the bidimensional nature of the system, the solid fraction φ is defined as the area fraction in the plane of shear, that is φ = Nπā 2 p /(L x h), withā p the average radius of a particle. Particles interact through the pairwise lubrication forces and frictional contact forces, with friction coefficient µ p = 0.5 (details of the interaction model can be found in [21]). Finally, particles are also coupled to the fluid via the Stokes drag. The fluid itself is unaffected by the particles, and follows a purely horizontal velocity field with x componentγy. The dynamics is overdamped (Reynolds and Stokes numbers are in the vanishing limit), so that the equation of motion corresponds to mechanical equilibrium on each particle in the bulk [22], which is solved for the velocities of the particles.
While the lower wall is fixed and not permeable to the fluid, the upper wall moves horizontally with a velocity u x =γh and is also free to move vertically with a velocity u y , according to the value of an imposed external vertical pressure P ext . The vertical motion implies a fluid flow through the permeable top wall, which imposes a Stokes drag on the particles making the wall. The vertical velocity of the wall is then determined at every time step by requiring that the upper wall is under mechanical balance, that is, imposing that the sum of the force applied by the bulk particles, the external force, and the Stokes drag from the background fluid (which is proportional to u y ) vanishes.
This setup mimics a constant external normal stress rheometer [4]. Just like in experiments, the bulk particle pressure P p equals the external one P ext (and accordingly J = J ext ) only in steady state. During transients however, J does not match J ext and moreover is not uniform in the sample, due to dilation effects as the particle pressure arises from both the imposed external pressure P ext but also from the fluid drag. Note that we do not take into account the fluid flow between pores which generates a drag on the particle phase during dilation or compaction. In our setup, dilation can only occur in the vertical direction, so that the problem can be reduced to a one-dimensional pore flow. This gives rise to a drag on the particle phase which is proportional to the difference of vertical velocities between fluid and particle phase, just like the Stokes drag, only with a smooth φ-dependent hydrodynamic resistance which is not singular at jamming. Hence, for the sake of simplicity our simulation framework ignores the φ-dependence of the hydrodynamic resistance of the pore flow, but does contain the same drag physics. This simplification does not affect the qualitative features of dilation, although at a quantitative level it underestimates the drag of the fluid phase on the particle phase during dilation, that is, our measured dilation dynamics is faster than it would be with full pore flow drag.

Transients following a step in J ext
We follow a systematic protocol to study the transients using a step in J ext , as illustrated in Fig. 3. The system is first brought to a steady state under an initial viscous number J i (black lines in Fig. 3), and at a strain γ = γ s = 50, J ext is suddenly changed to a value J f (red lines). It can be achieved in two ways: either by keepingγ constant while applying a step-change of P ext , or by keeping P ext constant while applying a step-change inγ, and here we adopted the former protocol. By dimensional analysis, both the protocols are strictly equivalent. If J f > J i , as in Fig. 3, the system dilates, as illustrated by the increasing height h as a function of the applied strain after the step in J ext . This dilation is associated to an instantaneous drop (resp. increase) of the fluid (resp. particle) pressure in the system, as well as a large pressure gradient, as seen in the difference between the particle normal stress P p normalized by applied pressure P ext measured at the bottom and top walls in Fig. 3d. Close to the bottom wall, the amplitude of the stress is reaching values two orders of magnitude above its steady state value. This value will depend on the sample height, however. Assuming a uniform dilation (which is actually not quite the case, see below), and using the Sus-  pension Balance Model [23][24][25], we expect the amplitude of the stress overshoot at the bottom of the system to increase roughly quadratically with the system height. The resulting fluid pressure gradient forces the fluid to permeate down through the particle phase, which in turn dilates. When dilation ends, P p /P ext reaches its steady state value of 1.
Dilation also has a strong effect on the macroscopic friction coefficient measured on the top wall µ wall = τ wall /P ext , with τ wall the shear stress on the wall, as shown in 3e. When J ext is increased, µ wall discontinuously increases by roughly two orders of magnitude before reaching its new steady state value µ SS (J f ). Hence the suspen-  sion transiently appears much more "frictional" than in steady state, due to dilation. Interestingly, dilation is not uniform across the height of the sample, as is evidenced by the snapshots of the system at different strains during dilation, as shown in Fig. 3a. Dilation occurs initially faster close to the top wall than in the bulk. This non uniform dilation can be related to the nonlinearity of the particle normal stress gradient, through the Suspension Balance Model [23][24][25]: at γ s , when the system is uniform in φ, SBM predicts that ∂ t φ(y) ∝ −∂ 2 y P p , so that the dilation is uniform if ∂ y P p is an affine function of y. As shown in the inset of Fig. 3d, right after the change of J ext , we see that this condition is not met. The stress gradient is superlinear close to the top wall, which leads to a faster dilation in this region.
Varying J f , as in Fig. 4(top), we see that dilation takes place on a strain scale that depends strongly on J f . Simulations with larger J f values take a much larger strain to reach steady state than simulations with a smaller J f . This is not only due to the fact that at larger J f values the system dilates by a larger amount, but is also a consequence of the smaller normal stresses generated by the suspension when its solid fraction decreases, which in turn lead to lower upper wall vertical velocities.
Comparatively, the effect of varying J i is much more short-lived, as shown also in Fig. 4(top): for a given J f , two simulations with J i = 1 × 10 −4 and J i = 1.6 × 10 −3 follow the same height evolution after a few strain units at most, often well before reaching steady state. The early transients however are very much affected by J i , and for a given value of J f one can observe dilation or compaction if, respectively, We model the dilation with a two-phase model, which consists of mass and momentum balances together with a constitutive model. Momentum balance, in the approximation of the Suspension Balance Model (SBM) [23][24][25], relates the particle normal stress along the velocity gradient σ p yy to the local particle phase velocity in the same direction u y as ∂ y σ p yy = 6η f φu y /ā p . ( Here we did not use the φ-dependent hydrodynamic resistance appearing in the SBM to be consistent with our simulation choice of a bare Stokes drag. In turn u y is linked to the dilation through mass conservation To close the problem, we use the steady µ(J) rheology [4] as a constitutive model, φ = φ SS (J), as fitted from the steady-state simulations (see Fig. 2). The boundary conditions require vanishing vertical velocity at the bottom wall u y (0) = 0 = ∂ y σ p yy (0), and force balance on the top wall, with u(h) = ∂ t h and N wall the number of particles in the upper wall. Altogether, this model gives the predictions shown in dashed line on Fig. 4(bottom), in good agreement with DEM results. Thus, in this setup extensions of the µ(J) rheology to capture transient aspects of the rheology [26] are not necessary to get a reasonable description of the dynamics. This is in contrast with the strong transients observed for the friction coefficient measured at the top wall µ wall : while globally the suspension is far from its steady state, a model assuming a local steady rheology is able to capture the dilation accurately. A possible rationalization of this result is that here dilation (or compaction) occurs on strain scales typically larger than 1, whereas the typical strain scale for a system to reach the steady µ(J) rheology is smaller than 1 (Pailha and Pouliquen estimate is ≈ 0.24 [26]). As a consequence for practical purpose the suspension can be considered locally in quasi steady state at all time.

Conclusion
We presented simulations of normal-stress controlled shear of granular suspensions, focusing on the transients during dilation. In particular, we capture the transient increase of macroscopic friction coefficient measured on the top wall. While the steady-state rheology is, as expected, well captured by the µ(J) rheology, the transient dynamics is also well predicted by assuming that the µ(J) rheology is locally followed at all time. We suggest that this will be true whenever dilation takes place over strain scales larger than 1, which would give the µ(J) rheology a vast domain of validity even in unsteady flows. Conversely, more stringent tests of the transient rheology, possibly showing limitations of the steady µ(J) and the need to extend it, should share the characteristics of generating quick volume changes. We are currently investigating this possibility with shear protocols differing from step changes in J ext , and we will report these results in a future article.