Viscometric flow of dense granular materials under controlled pressure and shear stress

This study examines the flow of dense granular materials under external shear stress and pressure using discrete element method simulations. In this method, the material is allowed to strain along all periodic directions and adapt its solid volume fraction in response to an imbalance between the internal state of stress and the external applied stress. By systematically varying the external shear stress and pressure, the steady rheological response is simulated for: (1) rate-independent quasi-static flow, and (2) rate-dependent inertial flow. The simulated flow is viscometric with non-negligible first and second normal stress differences. While both normal stress differences are negative in inertial flows, the first normal stress difference switches from negative to slightly positive, and second normal stress difference tends to zero in quasi-static flows. The first normal stress difference emerges from a lack of co-axiality between a second-rank contact fabric tensor and strain rate tensor in the flow plane, while the second normal stress difference is linked to an excess of contacts in the shear plane compared to the vorticity direction. A general rheological model of second order (in terms of strain rate tensor) is proposed to describe the two types of flow, and the model is calibrated for various values of interparticle friction from simulations on nearly mono-disperse spheres. The model incorporates normal stress differences in both regimes of flow and provides a complete viscometric description of steady dense granular flows.


Introduction
Granular flows exhibit several intriguing phenomena that distinguish them from Newtonian fluids, such as the presence of pressure-dependent arrest and flow onset (yield) criteria leading to rate-independent and rate-dependent flows, and a dilute gas-like flow dominated by inelastic particle collisions. A convenient classification defines three distinct types of granular flows (Forterre & Pouliquen 2008): (1) quasi-static flows, (2) dense inertial flows, and (3) gas-like collisional flows. The behavior of the three types of granular flows is quite diverse and several constitutive models have been proposed for their description; however, a general constitutive model applicable across all flow types remains elusive. In this work we focus our attention on the rate-independent and I. Srivastava, L. E. Silbert, G. S. Grest and J. B. Lechman rate-dependent flows, where the particle contact lifetimes are relatively long, inertia is important and the material is predominantly dense.
The rate-independent flow regime has been described by various constitutive models, largely inspired by the principles of solid mechanics and plasticity. Beginning with the incipient failure hypothesis of the Coulomb yield criterion (Sokolovskii 1965), further observations of critical state deformation in soils led to the development of several rigid-plastic and elasto-plastic models based on critical state theory and associated plasticity (Schofield & Wroth 1968). Recent developments have attempted to include material anisotropy in granular plasticity by introducing state-dependence of material stress, often characterized via material texture or fabric (Li & Dafalias 2012;Gao et al. 2014;Sun & Sundaresan 2011). The rate-independent granular plasticity has also been characterized by double shearing models (Spencer 1964) that relax the assumption of homogeneous deformations to explain shear banding along slip planes in granular materials. Further advances on these models have introduced the concepts of dilatancy (Mehrabadi & Cowin 1978), work hardening (Anand & Gu 2000), and more recently fabric anisotropy (Nemat-Nasser 2000; Zhu et al. 2006). The reader is referred to a recent review of various constitutive models of rate-independent regime in granular flows .
Rate-dependent granular flows were first observed to exhibit quadratic scaling of shear and normal stress with strain rate at a constant volume (Bagnold 1954), which was later verified in several experiments and simulations (Da Cruz et al. 2005;Lois et al. 2005;Pouliquen 1999;Silbert et al. 2001). Recently, a Bingham-type µ(I) rheological model for granular materials at moderate shear rates has been proposed (Jop et al. 2006), which attempts to connect the rate-independent and rate-dependent granular flow regimes by introducing pressure as a control variable instead of volume, although the two can be interchanged based on the one-to-one relationship between µ and I at moderate shear rates. Here µ is the dimensionless shear stress ratio, and I is the dimensionless inertial number (described later in the text). A detailed discussion of such viscoplastic models can be found in a recent review (Goddard 2014).
Although these rheological models have successfully predicted granular flow profiles in a remarkable number of geometries (MiDi 2004), several rheological effects remain unexplained, such as surface curvature in free-surface flows (Couturier et al. 2011;Mcelwaine et al. 2012), negative rod climbing effects (Boyer et al. 2011b), anomalous stress profile in Couette flows (Mehandia et al. 2012), and the observation of shear-free sheets in split-bottom Couette flows (Depken et al. 2007). Many of these effects arise from a lack of co-axiality between principal directions of stress and strain rate tensors in viscometric flows (Alam & Luding 2003, 2005Rycroft et al. 2009;Weinhart et al. 2013;Saha & Alam 2016;Bhateja & Khakhar 2018), which is the operating assumption in several constitutive models. As such, there is a need for higherorder constitutive models that incorporate these effects to provide better predictions of granular rheology. Furthermore, microstructural origins of these rheological effects, especially in the dense and quasi-static flow regimes, remain unclear.
In this paper, we describe fully stress-controlled discrete element method (DEM) simulations in both rate-independent and rate-dependent regimes. This simulation method enables the evolution of all strain degrees of freedom of a fully-periodic representative volume element of granular material in response to external applied shear stress and pressure. The novelty of this simulation method is four-fold: (1) it naturally captures the pressure-dependent flow-onset (yield) and flow-arrest phenomena (Srivastava et al. 2019b), (2) by prescribing shear stress rather than shear rate, this method can seamlessly traverse across rate-dependent and rate-independent regimes of granular flow, (3) by prescribing pressure rather than solid volume fraction, shear-induced dilation of granular materials is fully captured, and (4) the fully periodic nature of the simulations is devoid of any boundary effects, and thus represents the true bulk response of granular materials to applied stresses. The stress-controlled method is used to simulate shear flows of nearly mono-disperse spheres. A second-order rheological model that does not assume co-axiality of stress and strain rate tensors is proposed, and the model is calibrated for various values of interparticle friction from the simulation data. For viscometric flows, the second-order rheological effects result in non-negligible first and normal stress differences. We provide microstructural insights into the origin of normal stress differences by analyzing a secondrank contact fabric tensor.
The paper is organized as follows. Section 2 on Model and Methods describes the second-order rheological model, introduces the stress-controlled DEM simulation method, and provides details on the sphere-sphere contact mechanics model and general simulation parameters. Section 3 provides evidence that the steady flow generated in these simulations by applying shear stress and pressure is viscometric in nature. Section 4 describes the calibration of the rheological model based on the viscometric flow data from DEM simulations. Section 5 describes the normal stress differences measured in these simulations and their microstructural origins.

Rheological Model
In anticipation of the results presented below, we introduce a purely-dissipative rheological framework proposed by Goddard (1984), which was utilized to formulate constitutive laws for rate-independent and rate-dependent flow in granular materials (Goddard 1986(Goddard , 2014. In this framework, stress emerges from dissipation through macroscopic bulk deformation, which dominates over grain-scale inertial relaxation. This is similar to the dense flow of granular materials at low inertial numbers, which is of interest here. Furthermore, elastic effects are ignored, and non-hydrostatic stress components emerge entirely from granular flow, and are zero when there is no flow. In this framework, the Cauchy stress tensor σ is given by: where D is the symmetric strain rate tensor, and η{H} is a positive-definite fourth-rank tensor adhering to the constraints of a purely dissipative material, i.e., D : η{H} : D > 0, and is dependent upon the history H of flow, which can be conveniently represented through a deformation gradient F relative to a reference state. For the specific case of |D| → 0 corresponding to rate-independent plastic deformation (here, |D| = 1 2 D : D), the stress is given as: where µ 0 {H} is a fourth-rank yield modulus. Therefore, the total stress can be partitioned into its rate-independent and rate-dependent components as: where η 0 {H} is a fourth-rank viscosity tensor. We adapt this rheological framework to model steady granular rheology in our simulations through the following simplifications: (1) the steady values of η 0 and µ 0 are assumed to be independent of history H, and (2) η 0 and µ 0 are assumed to depend on D and its invariants. This dependence is introduced in a frame-indifferent manner to produce a second-order rheological model that well-describes non-isotropic flow effects observed in our simulations. With these simplifications, σ can be represented as: where p = 1 3 tr (σ) is the isotropic pressure. The deviatoric stress, σ − pI, depends on D and a spin tensor W . Here,γ = |D| is the magnitude of the strain-rate tensor, and |W | = 1 2 W : W is the magnitude of the spin tensor. The second, third and fourth terms in (2.4) represent rate-dependent contributions to the total stress that are characterized by the flow functions η 1 (γ, p) and η 2 (γ, p) and η 3 (γ, p), and are similar in form to a secondorder description of non-Newtonian fluids using Rivlin-Erickson tensors (Rivlin 1955). The fifth, sixth and seventh terms in (2.4) represent rate-independent contributions to the total stress that are characterized by plastic yield-like functions κ 1 (p), κ 2 (p) and κ 3 (p) that generally depend on the flow history. The pressure dependence of the flow functions is similar in spirit to the implicit constitutive theory of Rajagopal (2006). In this work we focus on shear flows, but in general, the coefficients η 1 , η 2 and η 3 also depend on tr(D 3 ), which is important when modeling non-viscometric flow scenarios . Similarly, the coefficients κ 1 , κ 2 and κ 3 depend on tr(D 3 ) |D| 3 , and this dependence can be derived from anisotropic models of granular plasticity. The contributions to the total stress from η 3 and κ 3 originate from the rotation of the principal directions of D with respect to σ, and the corresponding terms in (2.4) represent the steady state components In this paper, we will consider simple steady shearing flow of granular materials resulting from constant applied external shear stress and pressure, in which the memory of the flow has decayed and the deformation history is unimportant, thereby considerably simplifying the rheological model. In simple shear, the bulk velocity gradient ∇u takes the following viscometric form: for a flow along y direction, velocity gradient along x direction, and vorticity along z direction. From the definition D = 1/2(∇u + ∇u T ), simple shear flow corresponds to tr(D 3 ) = 0. In such viscometric flows, η 1 , η 2 , and η 3 represent the standard viscometric flow functions for non-Newtonian fluids (Coleman et al. 1966) corresponding to shear stress, second normal stress difference and first normal difference respectively. Similarly, κ 1 , κ 2 and κ 3 represent the rate-independent analogues of the these flow functions.
Correspondingly, for such viscometric flows, the stress tensor takes the following general form:

6)
Viscometric flow of dense granular materials 5 where σ xx = σ yy = σ zz . Previous simulations on sheared granular flows have proposed a similar form for the stress tensor, such as in granular flows down an incline (Silbert et al. 2001;Weinhart et al. 2013), free surface flows (Mcelwaine et al. 2012), and in shearfree sheets model (Depken et al. 2006) that proposed σ xx = σ yy = σ zz for quasi-static granular flows in a split-bottom Couette cell (Depken et al. 2007). The rheological model reduces to the well-known µ(I) relationship (Jop et al. 2006) for sheared granular flows when the second-order coefficients η 2,3 = 0 and κ 2,3 = 0. In this case, the stress tensor is assumed to be co-axial with the strain rate tensor, and the two are related to each other by a scalar relationship: where the stress ratio µ = |σ − pI|/p and the inertial number I = |D|a/(p/ρ) 0.5 , for an average particle diameter a and material density ρ. The µ(I) function is related to the flow coefficients of the rheological model in (2.4) through:

Constant Stress Simulations
Steady sheared flows can be simulated by applying a constant strain rate or a constant stress on the granular material. Previous simulations on granular flows have imposed a constant strain rate either through a solid wall-driven flow (Da Cruz et al. 2005;Koval et al. 2009;Kamrin & Koval 2014;Salerno et al. 2018) or by shearing the periodic simulation domain (Campbell 2002(Campbell , 2005Otsuki & Hayakawa 2011;Sun & Sundaresan 2011;Srivastava et al. 2019b). Wall-driven granular flows often result in flow localization near the walls (Shojaaee et al. 2012b), which requires careful calibration of wall properties to extract the bulk rheological properties (Shojaaee et al. 2012a;Schuhmacher et al. 2017). While a constant strain rate at the periodic boundaries can produce a viscometric flow field without walls (Campbell 2002(Campbell , 2005Peyneau & Roux 2008;Sun & Sundaresan 2011), it often results in large shear stress fluctuations (Peyneau & Roux 2008), especially in the quasi-static flow regime, which makes it challenging to calibrate the rateindependent part of granular rheology. Additionally, it was recently demonstrated that near the critical yield stress, granular flows are highly intermittent with a stochastic flow-arrest transition behavior. As such, a constant stress boundary condition is able to provide an accurate prediction of the rheology near the yield stress (Srivastava et al. 2019b). In this work, we simulate granular flows by applying a constant shear stress at the periodic boundaries, in which material is allowed to flow or not depending on the magnitude of applied stress. We will demonstrate that this boundary condition results in a well-defined viscometric flow.
Granular flows can also be simulated either at constant volume (isochoric) (Campbell 2002;Sun & Sundaresan 2011;Otsuki & Hayakawa 2011) or by imposing a constant normal stress (Campbell 2005;Sun & Sundaresan 2011;de Coulomb et al. 2017;Srivastava et al. 2019b). Granular materials dilate upon shearing, resulting in significant differences in the rheology between the two conditions (Campbell 2005). When the applied normal stress is constant, the material can dilate or compact upon shearing (depending on the initial condition) towards a 'critical state' solid volume fraction in the quasi-static regime (Schofield & Wroth 1968;Srivastava et al. 2019b). Furthermore, granular materials exhibit shear-induced dilation in the inertial regime. Isochoric granular flows are not commonly observed in practice, and various experiments often 6 I. Srivastava, L. E. Silbert, G. S. Grest and J. B. Lechman naturally correspond to a constant normal stress condition, such as in free surface flows (Mcelwaine et al. 2012;Jop et al. 2006) or flows in Couette cells (Lu et al. 2007;Dijksman et al. 2011). In this work, we simulate granular flows at a constant applied pressure where the material is free to adapt its volume. A constant pressure condition is different from case where all the normal stress components are specified equal to each other (Peyneau & Roux 2008).
To simulate the evolution of a granular system under constant external stress and pressure, we utilize a modularly-invariant adaptation (Shinoda et al. 2004) of the Parrinello-Rahman method (Parrinello & Rahman 1981) for molecular dynamics. This method was originally introduced to simulate the bulk properties of molecular systems in an isoenthalpic-isotension ensemble, including any phase transitions induced by the applied external stress (Parrinello & Rahman 1981). In the present simulations, a collection of particles contained within a 3D triclinic periodic cell is allowed to evolve under the application of a constant external stress tensor σ ext , which is constrained by (i) (1/3) σ ext,ii = p ext , (ii) σ ext,ij = τ ext for i, j = 1, 2 and 2, 1, and (iii) σ ext,ij = 0 for all other Einstein indices i = j. . Because the traction at the boundaries of the periodic cell is prescribed, the periodic cell itself is allowed to dilate (or compact), distort and rotate in all possible ways, thus simulating the true bulk response of the granular material under external stress. The triclinic periodic cell is represented by a cell matrix H which is a concatenation of the three lattice cell vectors that define the periodicity of the system. Upon the application of σ ext , the equations of motion for N particle positions and momenta {r i , p i }, and the periodic cell matrix and its associated momentum tensor {H, P g } are given by: where F i is the net force on a particle i, V is the variable volume of the periodic cell, and I is the identity tensor. A 'fictitious' mass W g associated with the inertia of the periodic cell is set as W g = N k n a 2 /ω 2 g , where a is the mean diameter of a particle, and ω g is the oscillation frequency associated with periodic cell fluctuations. The choice of ω g controls the magnitude of strain rate fluctuations during steady granular flow, but it does not affect its long-time mean value which is used to extract the rheology. We set ω g = 0.1ω p , where ω p = m e /k n is the frequency associated with the mechanical contact between two particles, and this choice provides a reasonable compromise between reducing the amplitude of strain rate fluctuations and reducing the simulation time required to attain steady flow under constant external stress.
The first two terms of the right side of (2.9d ) respectively represent the imbalance between bulk internal stress of the granular system σ int and external applied stress, which drives the motion of the periodic cell. The components of the bulk internal stress Viscometric flow of dense granular materials 7 σ int are calculated as: where r ij and F ij are the branch vector and the force between two contacting particles i and j, and v ′ i is the fluctuating velocity of particle i. Hereafter, the subscript 'int' will be dropped while referring to the internal state of the stress of the granular system. In (2.9d ), the tensor Σ is defined as: where H 0 is some reference state of the periodic cell, and J −1 HΣH T represents the 'true' measure of the external deviatoric stress, which is defined with respect to the reference state. Here J = det [F ] is the Jacobian of the deformation gradient F , which is defined as: As a result, this implementation of a constant external stress on the granular system prescribes the second Piola-Kirchoff measure of the external stress, or equivalently the thermodynamic tension (Souza & Martins 1997). In the present simulations, the reference state is updated to the current state at the end of every time step of integration of the equations of motion, in order to minimize the deviation of internal strain energy from work done by the external stress. All the simulations are performed using the large-scale molecular dynamics software LAMMPS (Plimpton 1995). The stress-controlled protocol described above has been previously implemented to study jamming (Dagois-Bohy et al. 2012;Smith et al. 2014) and creep (Srivastava & Fisher 2017) in granular packings, and recently to analyze flow-arrest transition in granular flows (Srivastava et al. 2019b).

Bulk Deformation
Upon applying an external pressure p ext and shear stress τ ext to a granular system, all the components of the macroscopic internal stress tensor σ evolve independently with time. Correspondingly, the triclinic periodic cell-represented by the matrix H-also evolves with time, and any changes in H can be decomposed into bulk volumetric and shear deformation, and rigid-body rotation. Consider the position r(t) of a particle at a simulation time t within the periodic cell, defined with respect to an origin (typically, one of the corners of the periodic cell). Its reduced coordinates s(t) can be defined by: (2.13) such that 0 < s(t) < 1. The periodic tiling of the space by the triclinic cell H(t) implies that a spatial coordinate r ′ (t) = H(t)[s(t)+∆] represents the periodic image of r, where ∆ is a vector of integers. The velocity v(t) =ṙ(t) of the particle is defined such that: (2.14) where the first term represents the contribution from the bulk periodic cell deformation and the second term represents the fluctuating non-affine velocity. Consequently, a bulk velocity gradient can be defined as ∇u(t) = ∇ r Ḣ (t)s(t) . Upon substituting (2.13): (2.15) 8 I. Srivastava, L. E. Silbert, G. S. Grest and J. B. Lechman The bulk velocity gradient can be further decomposed into a symmetric strain rate tensor D(t) = 1/2 ∇u(t) + ∇u T (t) representing the rate of volumetric and shear deformation, and an anti-symmetric spin tensor W (t) = 1/2 ∇u(t) − ∇u(t) T representing the rate of rigid-body rotation.

Contact Mechanics
In the present simulations, frictional spherical particles interact only upon contact. The contact forces are modeled using a spring and a dashpot along with a static yield criterion to model contact friction. This model was first developed by Cundall & Strack (1979), and since has been tested and implemented in various granular flow simulations (Silbert et al. 2001;Campbell 2005;Rycroft et al. 2009;Sun & Sundaresan 2011).
is the vector connecting their centroids. The contact normal force F nij and tangential force F tij on particle i are given by: where k n,t and γ n,t are elastic and damping constants, and m e = m i m j /(m i + m j ) is the effective mass. The corresponding force on particle j is given by Newton's third law such that F ji = F ij . The unit normal along the axis of contact is given by n ij = r ij /|r ij |, and v nij and v tij are respectively the normal and tangential components of the relative An elastic displacement u tij representing shear in the tangential direction is tracked during the lifetime of a contact, and it evolves according to the following ODE: with u tij = 0 at the initiation of the contact. Tangential friction between two contacting particles is modeled by a static yield criterion |F nij | < µ s |F tij |, which is always satisfied by limiting the tangential shear displacement u tij . The particle coefficient of friction µ s is a measure of its surface roughness, and significantly impacts the rheology of granular flow. The normal and tangential viscous damping at a contact are controlled by the coefficients of restitution e n,t = exp(−γ n,t t c /2), where t c = π(2k n /m e − γ 2 n /4) −1/2 is the characteristic time of contact between two particles (Silbert et al. 2001).

Simulation Details
The contact stiffness between particles k n and k t are set to unity. The normal damping constant γ n = 0.5 and the tangential damping constant γ t = 0.5γ n . Initially, dilute configurations of granular systems at a solid volume fraction φ = 0.05 are subjected to a constant external shear stress and hydrostatic pressure. We simulate granular flow at three external pressures p ext a/k n = 10 −4 , 10 −5 , 10 −6 , all within the asymptotic limit of rigid particle regime ( from τ ext /p ext = 0.0 to τ ext /p ext = 1.2 to simulate flows at various shear rates, and three different realizations are simulated for each shear rate. Each simulation consists of N ∼ 10 4 frictional particles whose diameters are uniformly distributed between 0.9a and 1.1a. Several particle coefficients of friction ranging from µ s = 0.0 to µ s = 1.0 are analyzed to study the effect of friction on stress-controlled granular rheology. Contact mechanics between two particle is resolved by setting the simulation time step to 0.02t c , where t c is the smallest characteristic timescale in the simulation. In the results presented below, time is scaled by t c , length is scaled by a, energy is scaled by k n a 2 , and stress is scaled by k n /a.

Evolution Towards Viscometric Flow
When the external shear stress τ ext and pressure p ext are switched on at t = 0, a dilute assembly of particles at an initial solid volume fraction φ = 0.05 responds with rapid volumetric compaction, as shown by the evolution of φ in figure 1(a) for a particular case of interparticle friction µ s = 0.3, p ext = 10 −4 and τ ext = 5 × 10 −5 . To estimate the total deformation accumulated by the material beyond isotropic compaction, we calculate the deformation gradient F (t) = H(t)H −1 0 as defined in (2.12), where H 0 is the periodic cell at t = 0. The rapid volumetric compaction at early times is seen by an equivalent decrease in F ii in figure 1(b) for i = x, y, z. The shear component F xy exhibits a super-linear increase at early times as a result shear deformation at low solid volume fractions in the absence of any significant resistance to the applied shear. At long times, F xy increases linearly with time, while F ii is constant and the other two shear components are negligible, thus indicating the achievement of steady incompressible viscometric flow. Although our focus here is on steady flows, the simulation method and rheological analysis described above provide the capability to calibrate a general historydependent rheological model defined in (2.3) for transient granular flows under constant or time-varying applied stresses.
Further insight into the nature of viscometric flow is given by the eigenvalue decomposition of the symmetric tensor D(t). We use the convention that the three orthonormal eigenvectors of D:d 1 ,d 2 andd 3 are ordered in the decreasing order of signed eigenvalues d 1 , d 2 and d 3 . Figure 1(c) shows the evolution of d 2 /d 1 and d 3 /d 1 as a function of time. At early times, the sum of eigenvalues is positive, which corresponds with rapid volumetric compaction as described above. The long time steady state flow is characterized by d 3 = −d 1 and d 2 = 0, which is a signature of planar flow, where the flow plane is spanned byd 1 andd 3 . To further ascertain the nature of planar flow, we calculate a vorticity parameter β defined as: where G =d 3d1 −d 1d3 . Figure 1(d ) shows the evolution of β with time. When the system transition into steady state flow at long times, β = 1, indicating simple shear deformation in the flow plane, thus confirming the viscometric nature of flow. During the transient evolution at early times, 0 < β < 1, indicating a complex flow behavior that is a mix of vorticity-free elongational flow (β = 0) and simple shear flow (β = 1) (Wagner & Mckinley 2016;). However, the flow at all times is homogeneous within the periodic cell and not localized. In steady state, the bulk rheological quantities fluctuate around their mean values, as seen in figure 1(a-d ). In order to achieve robust statistics, every simulation is run for at least 10 7 time steps to guarantee the achievement of steady state flow. This is especially necessary near the critical yield stress, where steady state equilibration takes a long time. Upon achieving steady state, each simulation is continued to run for at least another 10 6 time steps, during which the all data of interest are recorded at every 10 time steps and averaged to estimate their steady mean value. The statistical uncertainty associated with mean estimation is measured by its standard error using a block averaging method (Flyvbjerg & Petersen 1989). This method not only provides robust estimates of uncertainty around a mean value, but also indicates if the data has any long-time correlations, which would necessitate longer simulation runs for meaningful averaging.

Model Calibration
In this section, friction-dependent functional forms of all flow coefficients in (2.4) will be described. The material constants associated with these flow coefficients are extracted from the DEM simulation data by utilizing the fact that the four tensors The two flow coefficients η 1 and κ 1 have a first-order contribution (in terms of D) to the total stress σ, and they provide a measure of the shear stress in viscometric flow. These coefficients are estimated by: where τ = 1 2γ σ : D is the total flow-induced shear stress in the system.
Previous research has shown that shear flow of granular materials can be well-described by a local rheological relationship between a stress ratio µ and an inertial number I (Jop et al. 2006). In the present model, the stress ratio (hereby written with a subscript 1) is µ 1 = (η 1γ +κ 1 )/p, where η 1γ /p is the rate-dependent contribution and κ 1 /p is the rate-independent contribution. As such, η 1 represents the effective shear viscosity and κ 1 /p represents the static yield coefficient. Figure 2(a) shows the variation of µ 1 with I for five interparticle friction µ s at three p ext . All the curves at various pressures collapse onto a master curve for each µ s , which can be approximated by the following functional form described previously by Kamrin & Koval (2014): where, µ 0 1 , µ 1 1 , A and B are fitting parameters (see the table in the Appendix for the friction-dependent values of A and B). In the quasi-static, rate-independent regime wherė γ → 0, the stress ratio reaches a constant value µ 1 → µ 0 1 , which is equivalent to κ 1 /p in the rheological model. Although the rheology is unaffected by the applied pressure, as also observed in (de Coulomb et al. 2017), lower values of inertial numbers are achieved when the confining pressure is low, as seen by the variation of µ 1 − µ 0 1 with I for three pressures and two µ s in figure 2(b). Previous pressure and shear rate controlled simulations had demonstrated that the transition from geometric to kinetic flow regimes occurs at lower inertial number for lower confining pressure (de Coulomb et al. 2017), thus confirming the current observations. However, the present simulations produce highly stochastic flows in the quasi-static regime, which often arrest in the vicinity of the static yield coefficient (Srivastava et al. 2019b). Therefore, the rheology at low inertial numbers is not well-resolved for low pressures, especially for intermediate interparticle friction, as seen in figure 2(a). Recent experiments and simulations have indicated that the local rheology of frictional granular materials possibly exhibits hysteresis at very low inertial numbers (Degiuli & Wyart 2017), which would prohibit very slow flows in the present stress-controlled simulations.
The quasi-static stress ratio increases with friction from µ 0 1 = 0.13 for frictionless particles to µ 0 1 = 0.36 for particles with high friction, as shown in figure 2(c). The non-zero value of µ 0 1 for frictionless particles is consistent with previous simulations that demonstrated a non-zero internal friction angle for frictionless granular system (Peyneau & Roux 2008). The value of µ 0 1 at large friction is consistent with previous simulations and experiments (Boyer et al. 2011a;Salerno et al. 2018;Srivastava et al. 2019b), and is also similar to the critical stress ratio from the critical state theory (Schofield & Wroth 1968). For frictionless particles, µ 0 1 is slightly larger than µ 0 1 ∼ 0.11 reported in previous simulations (Peyneau & Roux 2008). This discrepancy could be attributed to the fitting errors arising from a lack of available data at low inertial number in the present simulations. The variation of µ 1 1 with µ s is similar to variation of µ 0 1 , and µ 1 1 > µ 0 1 , indicating a transition from the linear relationship between µ 1 and I at high inertial numbers to the quasi-static regime at low inertial numbers.
The flow coefficient η 1 , which corresponds to a shear viscosity, can be re-scaled to its dimensionless value as p −1/2 η1 d √ ρ . Figure 2(d ) shows the two friction-dependent limits of viscosity η 0 1 and η ∞ 1 corresponding to quasi-static regime and high shear rate flows respectively. The low shear rate viscosity η 0 1 decreases with increasing µ s , whereas the high shear rate viscosity η ∞ 1 is independent of µ s . At low shear rates, the average lifetime of a contact between two particles is long enough that friction crucially dominates the rheology of a granular system. However, at high shear rates, the average lifetime of a contact is short, and the effective shear viscosity is independent of friction.

4.2.
Flow Functions: η 2 and κ 2 In addition to the shear stress contribution to the total internal stress, there are nonnegligible second-order contributions that are typically observed in the flow of non-Newtonian fluids. In a viscometric description of such fluids, these effects are characterized by normal stress difference functions (Guazzelli & Pouliquen 2018). In the present rheological model, η 2 and κ 2 represent one set of such rate-dependent and rateindependent contributions. These coefficients are estimated by: and they represent the difference between mean normal stress in the flow plane and normal stress in the vorticity direction. A second stress ratio similar to µ 1 is defined as µ 2 = η 2γ 2 +κ 2 /p, where η 2γ 2 /p is the rate-dependent contribution and κ 2 /p is the rate-independent contribution. As such, η 2 represents the normal viscosity and κ 2 /p represents the normal yield coefficient. Figure 3(a) shows the variation of µ 2 with the square of inertial number I for five µ s and three p ext . All the curves at various pressures collapse onto a master curve for each µ s , which can be approximated by a functional form: (4.4) where, µ 0 2 , µ 1 2 , C and D are fitting parameters (see the table in the Appendix for frictiondependent values of C and D). In a manner similar to µ 1 , the quasi-static values of µ 2 at low inertial numbers are accessed for low confining pressures, as shown by the variation of µ 2 − µ 0 2 with I 2 in figure 3(b) for two different µ s that appear to collapse onto a single curve. However, the data at low inertial numbers is also noisy, resulting from the stochastic nature of slow granular flows.
In the quasi-static regime, µ 2 tends reaches a constant value µ 2 → µ 0 2 , which is equivalent to κ 2 /p in the rheological model. Its value varies monotonically from µ 0 2 = 0.03 for frictionless particles to µ 0 2 = 0.12 for particles with high friction, as shown in figure 3(c). The second fitting parameter µ 1 2 also increases monotonically with friction, and its value is always larger than µ 0 2 , indicating a transition from the linear relationship between µ 2 and I 2 at high inertial numbers to the quasi-static response at low inertial numbers.
Analogous to the shear viscosity η 1 , a normal viscosity η 2 is also extracted from the simulation data, which can be re-scaled to its dimensionless form: η 2 /d 2 ρ. Unlike η 1 , η 2 is independent of the applied pressure. Figure 3(d ) shows the two friction-dependent limits of dimensionless normal viscosity η 0 2 and η ∞ 2 corresponding to quasi-static and at high shear rate flows. The quasi-static viscosity is always larger in magnitude to the high rate viscosity for µ s . While η 0 2 is about three times larger for frictionless particles than particles with high friction, η ∞ 2 is largely independent of µ s .

4.3.
Flow Functions: η 3 and κ 3 An additional second-order contribution to the total stress emerges through the ratedependent and rate-independent flow coefficients η 3 and κ 3 respectively. They are estimated by: (4.5) and for viscometric flows, they represent the difference between the two normal stresses in the flow plane. A third stress ratio µ 3 is defined as: µ 3 = η 3γ 2 +κ 3 /p, where η 3γ 2 /p is the rate-dependent contribution and κ 3 /p is the rate-independent contribution. Here, η 3 represents an additional normal viscosity and κ 3 /p an additional normal yield coefficient. Figure 4(a) shows the variation of µ 3 with I 2 for five µ s at three p ext . All the curves at various pressures collapse onto a master curve for each µ s , which can be approximated by the following functional form: where µ 0 3 , µ 1 3 , E and F are fitting parameters (see the table in the Appendix for frictiondependent values of E and F ). The variation of µ 3 − µ 0 3 with I 2 exhibits similar trend with pressure in figure 4(b), as observed for µ 1 and µ 2 .
The stress ratio µ 3 exhibits a transition from negative values in the quasi-static regime to positive values at high inertial numbers for all µ s , as seen in figure 4(a) and also by the variation of µ 0 3 with µ s in figure 4(c). Unlike µ 0 1 and µ 0 2 , µ 0 3 is independent of µ s . The rate-dependent flow coefficient characterized by the normal viscosity η 3 can be re-scaled to its dimensionless form: η 3 /d 2 ρ. The two limits of dimensionless viscosity: η 0 3 in the quasi-static regime η 0 3 , and η ∞ 3 at high shear rates is displayed as a function of µ s in figure 4(d ). The quasi-static viscosity η 0 3 is approximately an order of magnitude larger than η ∞ 3 , and they both are about thrice as large for high friction particles compared to frictionless particles.

Granular Flow-Induced Dilation
The solid volume fraction φ of granular materials is highly sensitive to pressure and the rate of shear flow. These materials compact (jam) under the action of external pressure. However, under the action of external shear stress they dilate in order to flow, and the extent of dilation is higher for faster flows. In the present simulations, φ is not prescribed, and the system is allowed to freely attain its steady state solid volume fraction in response to the external stress. As such, we extract a dilatancy law relating the steady-state φ with the inertial number I of the flow. Figure 5 shows the variation of φ with I for five Figure 5. Solid volume fraction φ as a function of inertial number I for five interparticle friction µs (see legend) and three applied pressures: pext = 10 −4 , 10 −5 , 10 −6 . The vertical and horizontal error bars represent the standard error in the calculation of φ and I respectively. The black dashed lines represent fits to (4.7) for each µs. Inset: Variation of φ 0 and φ 1 with µs. µ s at three p ext . All the curves at various pressures collapse onto a master curve for each µ s , which can be approximated by the following functional form: φ 0 , φ 1 , G and H are fitting parameters (see the table in the Appendix for frictiondependent values of G and H).
The quasi-static solid volume fraction φ 0 varies significantly with µ s ranging from φ 0 = 0.637 for frictionless particles to φ 0 = 0.589 for particles with high friction, as shown in the inset of figure 5. This is consistent with the previous finding that frictionless particles do not dilate to flow (Peyneau & Roux 2008), and thus this value is similar to the random close packing of mono-disperse spheres. The φ 0 for high friction particles is also consistent with the critical state solid volume fraction from the critical state theory (Schofield & Wroth 1968). Additionally, the variation of φ 0 with friction is also consistent with previous experiments on granular flows (Boyer et al. 2011a) and simulations on flow-arrest transition in granular materials (Srivastava et al. 2019b).
The second fitting parameter φ 1 also varies similarly with µ s , and its value is always smaller then φ 0 , thereby indicating a transition from a constant quasi-static value to a linearly decreasing function with I at higher inertial numbers. The effect of µ s on steady φ decreases substantially at high I, and this can be attributed to short-lived contacts in fast flows that diminishes the effect of interparticle friction.

Normal Stress Differences and their Microstructural Origins
The non-negligible second-order contributions to stress in viscometric granular flows indicate the presence of normal stress differences. Previous research on sheared granular and suspension flows has demonstrated the existence of normal stress differences (Silbert et al. 2001;Alam & Luding 2005;Rycroft et al. 2009;Sun & Sundaresan 2011;Couturier et al. 2011;Boyer et al. 2011b;Weinhart et al. 2013;Saha & Alam 2016;Guazzelli & Pouliquen 2018), and these differences have been attributed to flow-induced microstructural effects in dilute granular flows (Saha & Alam 2016) and in dense suspension flows . Particularly, normal stress differences can arise either from: (1) a misalignment of σ and D in the flow plane, known as the first normal stress difference, or (2) from the anisotropy of normal stress between the flow plane and the vorticity direction, known as the second normal stress difference.
In this section, we describe the microstructural origin of normal stress differences in dense viscometric granular flows by utilizing a second-rank contact fabric tensor A, which provides a convenient description of the directional distribution of particle contact network and inherent structural anisotropy (Oda 1982;Kanatani 1984). The orientational distribution P (n) of contact normal unit vectors n can be expanded to the second order in Fourier series as (Rothenburg & Bathurst 1989;Bathurst & Rothenburg 1990): where A is trace-free and symmetric. For dense granular materials where internal stress σ is dominated by particle contacts, it can be expressed as the following integral in the orientational space Ω (Rothenburg & Bathurst 1989;Bathurst & Rothenburg 1990;Srivastava et al. 2019a): where l 0 and f i are the average center-to-center distance and the i-th component of the normal force between two contacting particles respectively. This representation of the stress tensor ensures that σ and A have the same structure.
5.1. Second Normal Stress Difference A significant source of normal stress anisotropy emerges from the difference between the mean normal stress in the flow plane and normal stress in the vorticity direction, and is represented by the viscometric flow function N 0 /τ = σ zz − σxx+σyy 2 /τ , where y is the flow direction, x is the flow gradient direction and z is the vorticity direction, and the stresses are defined negative in the compressive sense. In the present simulations, N 0 /|τ | is computed by: In figure 6(a), N 0 /τ is plotted as a function of the distance to quasi-static solid volume fraction φ − φ 0 for five µ s . The negative value of N 0 /τ implies that there is more stress in the flow plane than in the vorticity direction. In the present simulations, an imbalance between external applied pressure and internal pressure drives isochoric periodic cell deformation, as shown by the equal values of F ii (t) in figure 1(b). In the scenario where each σ ii is individually balanced, we observed a rapid compaction of the box in the vorticity direction leading to simulation instability arising from the second normal stress difference. The magnitude of second normal stress difference is larger for frictional particles than for frictionless particles; however, even frictionless particles exhibit non-zero second normal stress difference during flow. As the solid volume fraction increases towards quasi-static φ 0 , the anisotropy consistently decreases for all µ s . For frictionless particles, N 0 appears to tend to zero in the quasi-static limit corresponding to φ 0 = 0.637, which is the the random close packing density for mono-disperse spheres. The small non-zero value at the lowest inertial number could be a consequence of finite particle stiffness  or finite system size. However, the out of flow plane stress anisotropy is demonstrably non-zero for frictional particles in the quasi-static limit. The values of N 0 /τ in figure 6(a) match well with the recently reported values in simulations on inertia-less frictional suspensions . The notion of non-zero anisotropy in the quasi-static regime is also consistent with recent experiments on suspensions that predict the existence of a normal yield stress (De Cagny et al. 2019).
An implication of these findings is that the flow of frictional granular materials is not co-directional, i.e., the hypothesis σ ∝ D, which has been been assumed within the µ(I) rheological model is not accurate. N 0 /τ increases with I for all µ s , as seen in figure 6(b), and remains measurably non-zero for frictional particles even at low I. Prior simulations on quasi-static simple shear granular flows (Sun & Sundaresan 2011), granular flows down an incline (Silbert et al. 2001), gravity-driven granular flows through an orifice (Rycroft et al. 2009), and granular flows in a split-bottom Couette cell (Depken et al. 2007) have questioned the co-directionality hypothesis. Two previously proposed theoretical models-double shearing (Spencer 1964) and shear-free sheets (Depken et al. 2007)-have incorporated these effects for quasi-static and dense granular flows.
The second normal stress difference results from an excess of contacts oriented in the flow plane than those oriented in the vorticity direction. Figure 6(c) shows the variation of N a 0 /Z 2 with φ − φ 0 for various interparticle friction. Here, N a 0 = A zz − Axx+Ayy 2 is the contact fabric 'second normal difference', which represents the anisotropy in average orientation of contacts between the flow plane and the vorticity direction. The rattlerfree coordination number is computed as Z 2 = 2N c /(N − N r ), where N r is the number of rattler particles with less than two contacts, and N c is the total number of contacts with non-zero normal force belonging to non-rattler particles (Sun & Sundaresan 2011). At low φ corresponding to large I (see figure 6(d )), a higher fraction of contacts are oriented in the flow plane, which results in large normal stress difference. Upon approach to the quasi-static regime at high φ, the contact distribution becomes more isotropic, resulting in reduced normal stress difference. For frictionless particles, the orientational distribution of contacts expectedly becomes isotropic in the quasi-static regime, as seen by N a 0 → 0.

First Normal Stress Difference
The first normal stress difference, which characterizes the anisotropy between σ and D in the flow plane, is represented by the viscometric flow function N 1 /τ = (σ yy − σ xx ) /τ . In the present simulations, this is calculated from: The variation of N 1 /τ with φ − φ 0 for five µ s is displayed in figure 7(a). At low φ, N 1 is negative for all µ s , and its magnitude increases with increasing µ s for a given distance from the quasi-static solid volume fraction φ − φ 0 . At low φ, the ratio of N 0 /N 1 is approximately 3 − 4 for all interparticle friction, consistent with previous findings (Gallier et al. 2014).
When the flow becomes dense, N 1 increases towards zero and becomes positive for highly dense flows in the quasi-static regime. The change of sign of N 1 at high φ has been previously observed in simulations (Alam & Luding 2005;Weinhart et al. 2013;) and experiments (Couturier et al. 2011), but its existence is debated, and has been attributed to interparticle friction (Dbouk et al. 2013) and boundary wall effects in experiments (Gallier et al. 2014). In the present simulations, we observe slightly positive N 1 for all values of µ s at large φ, and there are no boundary effects in these bulk simulations. Recently it was demonstrated that finite particle stiffness-which is often used as numerical regularization in hard particle simulationscauses N 1 to become positive at large φ in simulations on inertia-less frictional suspensions . However, our granular simulations do not provide any conclusive evidence of vanishing positive N 1 as a result of increasing particle stiffness. A careful analysis about this effect constitutes an important part of our future work.
The first normal stress difference is related to the angular misalignment θ c between the principal directions of D (d 1 andd 3 ) and A (â 1 andâ 3 ) in the flow plane, as described in the schematic in figure 7(b). We setâ 1 andâ 3 in the compression and expansion direction of shear flow respectively. A negative θ c indicates that â 1 ×d 1 ·d 2 > 0.
The misalignment angle θ c and N 1 /τ are strongly correlated, as depicted in figure 7(c) for various µ s that largely collapse onto a single curve. The negative value of θ c results in excess stress along the gradient direction as compared to the flow direction, which sets the negative sign of first normal stress difference. Such microstructural origins of N 1 , arising from a misalignment between the projected contact vectors and principal flow direction in the flow plane, were previously demonstrated for inertia-less frictional suspensions , and in the case of dilute granular flows through a similar misalignment between fluctuating velocity moment tensor and principal flow direction in the flow plane (Saha & Alam 2016). Remarkably, θ c → 0 as N 1 → 0, indicating a vanishing misalignment between D and A at high solid volume fractions. This is also seen by the variation of θ c with I in figure 7(d ), where the data for all µ s collapse onto a single curve, and this variation is consistent with prior simulations on granular flows down an incline (Weinhart et al. 2013). A small positive first normal stress difference exists at high solid volume fractions in the vicinity of yield stress despite a near-complete alignment of D and A in the flow plane, thus indicating that either a different underlying physical phenomenon is responsible for positive N 1 , or it is possibly a consequence of finite system size.

Conclusions
In this paper, we described a discrete element method to simulate dense granular flows under external applied stress in a fully periodic representative volume element. Rather than prescribing solid volume fraction and/or strain rate, this method enables independent evolution of solid volume fraction and 3D strain rate tensor in response to an imbalance between internal state of stress and external applied stress. Using this method, bulk viscometric granular flows were simulated under external pressure and shear stress, which was devoid of any boundary effects, and thus closely represented the boundary conditions often found in practice.
We developed a second-order rheological model to relate the internal Cauchy stress σ with the strain rate tensor D for various interparticle friction. The model considers both rate-dependent and rate-independent contributions to the total stress, where the latter is often described using models of granular plasticity. The rheological model wellpredicts the µ(I) rheology of granular materials. Additionally, it also predicts normal stress differences in steady viscometric granular flows, which have often been observed in simulations and experiments, but have not been well-characterized. A major implication of this model is that it does not impose co-axiality between σ and D in dense granular flows, which is often assumed in several other constitutive models.
A major focus of this work has been to highlight the role of interparticle friction on viscometric granular rheology in the dense flowing regime, particularly on the two normal stress differences. We found that friction not only increases the quasi-static shear stress ratio, but also the quasi-static value of second normal stress difference, thus indicating the presence of an anisotropic yield stress. At higher flow rates in the inertial regime, friction consistently increases the magnitude of both normal stress differences, indicating an increasing departure from the co-axiality of σ and D. Although the second normal stress difference is always negative, the first normal stress difference changes its sign from negative to positive at high solid volume fractions in the quasi-static regime. Further microstructural investigations highlighted that negative first normal stress difference results from a misalignment between D and a second-rank contact fabric tensor A in the flow plane, which describes the orientational distribution of sphere-sphere contacts in granular flows. In contrast, the second normal stress difference results from an excess of contacts oriented in the flow plane than in the vorticity direction.
These results demonstrate the importance of developing rheological models beyond simple scalar models to predict granular rheology in complex flow fields, which are often observed in practice. The breakdown of co-axiality of stress and strain rate tensors highlights the need for a rheological model that includes contact fabric tensor as an internal variable. A general form of such models σ = F (D, A) was recently proposed for granular materials and suspensions (Goddard 2006;Stickel et al. 2006), and was calibrated for rate-independent granular flows (Parra & Kamrin 2019;Sun & Sundaresan 2011). The calibration of these models in rate-dependent flows is required. Additionally, extensions of these rheological models to non-viscometric flows such as uniaxial or triaxial compression, as well as transient evolving inertial granular flows is also required, and is currently a subject of our study.  Table 1. Fitting parameters corresponding to equations (4.2), (4.4), (4.6) and (4.7), as a function of interparticle friction µs.