In order to identify the parameters in the second-order model, we use LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator), which implements the discrete element method (DEM), to simulate many different homogeneous and inhomogeneous planar shear flows of 3D frictional spheres. The grain diameter d _i is set to be uniformly distributed from 0.8d to 1.2d to prevent crystallization. When we calculate I, the local average diameter is used. The mass of the grains is determined as m i = 4 3 π ( d i / 2 ) 3 ρ s , which distributes around m ≡ 4 3 π ( d / 2 ) 3 ρ s . For the contact forces, we use the standard spring-dashpot model introduced by Cundall and Strack [49] with Coulomb friction. This model has been used in many previous studies, such as to study simple shear flows in 2D [2] and various inhomogeneous flows in 2D [5, 17, 50, 51] and 3D [20, 25]. In this model, the normal force is given by F _n = k _n δ _n − γ _n v _n (repulsive) where k _n is the normal elastic constant, δ _n is the normal component of the contact displacement, γ _n is the damping coefficient, and v _n is the normal component of the relative velocity. We use k _n = 2.63 N/m × 10⁵ N/m, d = 0.0008 m, and ρ _s = 2500 kg/m³ which can be considered when reading the pressure data in Figure 13. The tangential force is given by F _t = −k _t δ _t for |k _t δ _t | < μ _p F _n and F _t = −μ _p F _n × (δ _t /|δ _t |) for |k _t δ _t | > μ _p F _n where μ _p = 0.4 is the surface friction coefficient and δ _t is the tangential component of the contact displacement obtained by integrating tangential relative velocity during the collision. The minus sign indicates the tangential force is pointing in the direction of decreasing δ _t . The tangential elastic constant k _t is set to be 2/7 of k _n such that the frequencies of normal and tangential contact oscillations are similar [43]. To simulate hard particles, the stiffness number k _n /Pd is kept larger than 10⁵. The damping coefficient for the restitution coefficient ϵ and the effective mass of collision m _eff, which is approximately m/2, is calculated as γ n = 4 m eff k n / ( 1 + ( π / ln ⁡ ϵ ) 2 ) [2, 51]. We set ϵ = 0.24 to match with Kim and Kamrin [25]. The choice of k _t and ϵ is known to have little impact on the flow behavior in the case of dense flows of hard particles, as reported in [2, 37, 52]. The simulation time step is set to be 6% of the binary collision time τ c = π m 2 k n 1 + ( ln ⁡ ϵ / π ) 2 which is half the period of an underdamped oscillator made of the two identical grains.

As in Kim and Kamrin [25], we test simple shear flows (Figure 2A), shear flows under gravity (Figure 2B), flows in vertical and tilted chutes (θ = 90° and 60° in Figure 2C), and “concave flows” (Figure 2D). The concave flow is named so because the plot of the shear rate against height becomes concave upward due to the outward body force F ⃗ z = ( m i G / d ) ( z − z 0 ) z ̂ for the midpoint of the system z ₀. G is a constant. The horizontal length and depth of the system are L _x = 20d and L _y = 16d respectively and the side boundaries are set to be periodic. Following previous studies [2, 5, 20, 25, 50, 51], we move the top wall in the z-direction to maintain the top-wall pressure P _wall. The horizontal wall velocity V _wall is set to be constant. More detailed DEM simulation conditions are in the Supplementary Material.

The continuum variables at each height are obtained by coarse graining, in the same way as [25]. The instantaneous packing fraction at time t can be calculated as ϕ z k , t = ∑ i A k i A (6) where A is the area of the horizontal plane and A _ki is the cross-sectional area between the ith particle and the plane of z = z _k . We kept the interval of z _k less than 0.5d. The macroscopic velocity field can be calculated as v ⃗ z k , t = ∑ i A k i v ⃗ i t ∑ i A k i (7) where v ⃗ i is the velocity of the ith particle. We define the instantaneous granular temperature tensor as T z k , t = ∑ i A k i δ v ⃗ i t ⊗ δ v ⃗ i t ∑ i A k i (8) where δ v ⃗ i ( t ) = v ⃗ i ( t ) − v ⃗ ( z k , t ) is the instantaneous velocity fluctuation of the ith particle. The instantaneous stress is σ z k , t = σ K z k , t + ∑ i A k i σ i t A (9) where σ _i is the contact stress of the ith particle, and σ ^K = −ρ _s ϕ T is the kinetic stress [43]. The specific formula for σ _i is σ i = 1 V i ∑ j ≠ i 1 2 r ⃗ i j ⊗ f ⃗ i j (10) where V _i is the volume of the ith particle, r ⃗ i j is the displacement from the center of the ith particle to the center of the jth particle, and f ⃗ i j is the interaction force exerted on the ith particle by the jth particle. For one particle, the stress tensor may not be symmetric. However, if a sufficient number of particles are averaged, the coarse-grained stress tensor becomes symmetric as shown by [53]. We choose ( T y y + T z z ) / 2 = ( δ v y 2 + δ v z 2 ) / 2 as the granular temperature T because for some flows, T _xx is measured quite differently from the other diagonal elements. The ratios between diagonal elements of the granular temperature in the inclined chute flow can be found in the Supplementary Material. These coarse-grained fields are averaged over time once the flow reaches a steady state. We cut off the data where total local shear is less than 1 or the distance from the walls is less than 3d to include enough configurations and exclude the wall effect.

The second-order model parameters μ ₁, μ ₂, and μ ₃ can be calculated from the measured stress tensor and velocity gradient tensor in the DEM planar shear flows (Figure 2). By symmetry, both homogeneous and inhomogeneous flows have negligible mean v _y and v _z and the mean v _x changes only in the z-direction. Therefore, the velocity gradient tensor has only one non-zero component: L x z = v x , z ≡ ∂ v x ∂ z which can vary depending on the height. With its magnitude γ ̇ , the inertial number can be obtained as I = γ ̇ / P / ρ s d 2 . By definition, the symmetric strain rate tensor D and the spin tensor W are D = 0 0 v x , z / 2 0 0 0 v x , z / 2 0 0 and W = 0 0 v x , z / 2 0 0 0 − v x , z / 2 0 0 (11) respectively. Using D 2 − t r D 2 3 I = v x , z 2 12 1 0 0 0 − 2 0 0 0 1 (12) and D W − W D = v x , z 2 2 − 1 0 0 0 0 0 0 0 1 , (13) we can express the stress tensor (Eq. 5) as σ = P − 1 + 1 3 μ 2 − 2 μ 3 0 μ 1 v x , z γ ̇ 0 − 1 − 2 3 μ 2 0 μ 1 v x , z γ ̇ 0 − 1 + 1 3 μ 2 + 2 μ 3 . (14) Because v x , z / γ ̇ = 1 for v _x,z > 0 and v x , z / γ ̇ = − 1 for v _x,z < 0, the sign of σ _xz is determined by the sign of v _x,z while the magnitude of σ _xz is μ ₁ P in either case. We can see that the second-order model parameters μ ₂ and μ ₃ introduce differences between the normal stress components. Inversely, the three parameters can be extracted from the stress measurement using the fact that the basis tensors I , D − 1 3 t r ( D ) I , D 2 − 1 3 t r ( D 2 ) I , and D W − W D are orthogonal to each other under the tensor inner product: μ 1 = 1 γ ̇ P σ : D = σ x z P v x , z γ ̇ = | σ x z | P μ 2 = 6 γ ̇ 2 P σ : D 2 − t r D 2 3 I = − σ x x + σ z z − 2 σ y y 2 P μ 3 = 1 2 γ ̇ 2 P σ : D W − W D = σ x x − σ z z 4 P (15) with the assumption that the flow follows the second-order fluid equation (Eq. 5).

The ratios between the coarse-grained normal stresses from the planar shear tests are plotted in Figure 3. Figure 3A shows that the normal stress in the flow direction σ _xx is nearly the same as the normal stress in the velocity gradient direction σ _zz for the low I regime (I < 0.1). As I increases beyond 0.1, the magnitude of σ _xx becomes larger than the magnitude of σ _zz and their ratio σ _xx /σ _zz keeps increasing. This behavior is similar to previous findings in 2D and 3D granular systems [35, 36, 38–40, 43, 45]. On the other hand, the normal stress in the vorticity direction σ _yy is spread between 85% and 95% of σ _zz for I < 0.1. As I increases beyond 0.1, the ratio becomes even smaller as shown in Figure 3B. This is in line with previous observations [37, 43, 45].

In the planar shear tests with one velocity component, the coefficient of the first-order derivative in the second-order model μ ₁ can be calculated as |σ _xz |/P (Eq. 15). Previously, Kim and Kamrin [25] have examined μ = |σ _xz /σ _zz | assuming P = −σ _zz and identified that μΘ ^p data collapse into a master curve f(I) when p ≈ 1/6. The second-order model, however, does not assume the normal stress isotropy and the DEM stress data actually exhibits significant anisotropy. Therefore, we measure μ ₁ = |σ _xz |/P with P = −(σ _xx + σ _yy + σ _zz )/3 and recalibrate the master curve that μ ₁Θ^1/6 data collapse into.

The scattered μ ₁ data in the μ ₁ vs. I plot in Figure 4A can be gathered to a single line by multiplying by Θ p 1 where p ₁ is roughly 1/6 as shown in Figure 4B, which is almost the same as the previous μ behavior in [25]: μ 1 Θ p 1 = f 1 I (16) where f ₁(I) can be fitted by a form I α 1 where α ₁ gradually increases from 0.25 to 0.5 as I increases from 10^–4 to 1. For a given I, f ₁(I) would increase with the surface friction μ _p as previous studies on homogeneous [45] and inhomogeneous [25] flows have suggested. For the finite difference method (FDM) simulations in Section 3, we use p ₁ = 1/6 and f ₁(I) = 0.141I ^0.21 + 0.132I ^0.4 + 0.29I ^0.8 − 0.050I ^1.6 which is arbitrarily chosen to depict the collapsed data of the μ _p = 0.4 case shown in Figure 4B.

The coefficient of the D 2 − t r ( D 2 ) 3 I term, μ ₂, can be calculated as −((σ _xx + σ _zz )/2 − σ _yy )/P according to Eq. 15. It represents a measure of the difference between the normal stress in the vorticity direction and the mean of the other normal stresses. Simply put, the larger the μ ₂, the smaller |σ _yy | is compared to the other normal stresses; Figure 5A shows that μ ₂ monotonically increases as I increases. This monotonic increase in the μ ₂ vs. I graph has also been observed by [45]. As I decreases below 10^–2, similar to μ ₁, the data points of μ ₂ in inhomogeneous flows become more scattered. The random error in μ ₂ measurement seems much larger than μ ₁ possibly because μ ₂ is measured from the small difference between (σ _xx + σ _zz )/2 and σ _yy accumulating the errors of the three stress elements. Although the μ ₂ data are not as clean as μ ₁, we can still see that μ ₂ exhibits Θ dependence similar to μ ₁ where higher Θ lowers μ ₂ for a given I. This may imply that μ ₂ can also be scaled by a power function of Θ like μ ₁. Indeed, multiplying μ ₂ by Θ p 2 seems to cancel the non-local spread and collapse the data more onto a single curve when p ₂ ≈ 1/6 as shown in Figure 5B: μ 2 Θ p 2 = f 2 I (17) where f ₂(I) can be fitted by I α 2 where α ₂ gradually increases from 0.3 to 0.7 as I increases from 10^–3 to 1. For a fixed I, f ₂(I) would increase with μ _p ; [45] has similarly shown for homogeneous flows that μ ₂ increases with μ _p . For the FDM simulations in Section 3, we use p ₂ = 1/6 and f ₂(I) = 0.093I ^0.29 + 0.195I ^1.2 − 0.035I ² (dashed line in Figure 5B) which is arbitrarily chosen to represent the collapsed data of the μ _p = 0.4 case.

As in [25], we plot μ ₂ vs. Θ at many fixed choices of I to see the functional dependence more clearly. We choose five panels of data in 0.9I* < I < 1.1I* for I* = {3 × 10^–4, 6 × 10^–4, 1.2 × 10^–3, 2.4 × 10^–3, 4.8 × 10^–3}. The dashed lines in Figure 5C illustrate the best-fit lines of a multivariate linear regression whose slope is measured as −0.154 ± 0.014 which equals −p ₂. The actual error could be larger than this standard error measured only from the chosen data. Since p ₂ seems not so different from 1/6 which is the exponent of μ ₁ and the data collapse is strong with 1/6 in Figure 5B, we assume both μ ₁ and μ ₂ scale with the same 1/6 power of Θ for the continuum simulations in Section 3. It remains for further research to identify the exponents and the master curves more accurately.

The coefficient of the DW − WD term, μ ₃, can be calculated as (σ _xx − σ _zz )/(4P) according to Eq. 15. It represents a measure of the difference between the normal stresses in the flow direction and the velocity gradient direction; Figure 6A shows that μ ₃ is near or slightly above zero for I < 0.1. The only exception is the data from the chute flow geometries, which go up to 0.02. This behavior is significantly different from the non-local effect of granular temperature observed in Figure 4A because the data from the other inhomogeneous flows (the shear with gravity and the concave flows) follow the same curve as the homogeneous flow data, which means flows with different Θ can have the same μ ₃ for a given I. Θ is not enough to explain this deviation and there must be other factors that affect the μ ₃ measurement in the chute flows. We discuss another possible μ ₃ calibration to resolve this issue in the next section. Here, we choose to exclude the problematic chute flow data in calibrating μ ₃. Figure 6B shows that, as I increases, μ ₃ becomes negative and decreases almost linearly: μ 3 = f 3 I (18) where f ₃(I) ≈ − 0.045I (dashed line) for μ _p = 0.4. For a fixed I, f ₃(I) would decrease (increase in magnitude) with μ _p as [45] has shown for homogeneous flows.

The chute flows’ peculiar μ ₃ behavior in Figure 6A may be attributed to the anisotropy of the granular temperature. This possibility suggests an alternative way to calibrate μ ₃ which could be implemented in future work. Let us denote the dimensionless granular temperature tensor as Θ i j = ρ s T i j P . (19) where T _ij is the (i, j)th element of the granular temperature tensor introduced in Eq. 8. The difference between normal elements in the flow direction and the vorticity direction, Θ _xx − Θ _yy , behaves similarly to μ ₃ in that the chute flow data diverge from the other data as I decreases as shown in Figure 7A. The anisotropy of the granular temperature tensor may cause the peculiar behavior of μ ₃, or there is another unknown macroscopic variable that affects both Θ _xx − Θ _yy and μ ₃ in a similar pattern. In any case, we can calibrate μ ₃ using this correlation. For example, adding an arbitrary function 0.4 ( Θ x x − Θ y y ) 1 / 6 to −μ ₃ seems to collapse the data more: − μ 3 + 0.4 Θ x x − Θ y y 1 / 6 ≈ f 3 A I (20) as shown in Figure 7B. This collapse only shows the correlation between μ ₃ and Θ _xx − Θ _yy in our DEM data. The actual (more accurate) expression for μ ₃ may differ from Eq. 20. More diverse Θ _xx − Θ _yy vs. I curves are needed to clearly verify the data collapse as our current flow geometries have generated only two different branches as can be seen in Figure 7A. Also, it is not obvious how to define Θ _xx and Θ _yy in a general flow. Therefore, we leave this problem for future research and choose to use the simpler μ ₃ expression, Eq. 18, for the continuum simulations in Section 3.

In this section, we examine the first and second normal stress differences as they are commonly measured quantities to represent the normal stress anisotropy of a material, even though we do not utilize them explicitly in our continuum model. The first normal stress difference is defined as N ₁ = σ _xx − σ _zz which is the same as 4μ ₃ P (Eq. 15). Figure 8A shows that N ₁/P is almost zero for I < 10^–1 and grows negative as I increases above 10^–1. It means that, as I increases, the magnitude of stress in the flow direction |σ _xx | becomes larger than that in the gradient direction |σ _zz |. This is consistent with previous observations in 2D and 3D granular systems [35, 36, 38–40, 43, 45]. The peculiar chute data do not collapse into a single line even if the horizontal axis is changed to |σ _xz /σ _zz | as can be seen in Figure 8B.

On the other hand, the second normal stress difference N ₂ = σ _zz − σ _yy represents normal stress anisotropy in the plane formed by the velocity gradient and the vorticity directions. N ₂ can be written as −(μ ₂ + 2μ ₃)P (Eq. 15). Figure 8C displays that N ₂/P behaves like μ ₁ and μ ₂ in that the inhomogeneous flow data are scattered and the simple shear data do not converge to zero in the quasi-static regime. However, the chute flow data are still a bit out of the trend. This peculiarity is more noticeable when the horizontal axis is |σ _xz /σ _zz | as shown in Figure 8D. All the other N ₂/P data seem to form a single line while the chute data with |σ _xz /σ _zz | < 0.4 have lower N ₂/P values. Another important feature of Figure 8C is that N ₂ is always negative for a non-zero I, which is in line with previous observations [37, 43, 45]. It has been known that a negative N ₂ makes the free surface convex up in a channel flow with no surface tension [29, 54–57]. This convex surface is also observed in our inclined chute flow simulations in Section 3.

If we look closely at Figures 8B, D where the horizontal axes are |σ _xz /σ _zz |, the chute flow data go higher than the other data in Figure 8B and lower in Figure 8D. Because the peculiar chute deviations have opposite signs in N ₁ and N ₂, we are intrigued to observe their sum N ₁ + N ₂. This variable is actually another normal stress difference σ _xx − σ _yy . Interestingly, (N ₁ + N ₂)/P exhibits a collapsed line when the horizontal axis is |σ _xz /σ _zz | canceling the above deviations as displayed in Figure 9B. Meanwhile; Figure 9A demonstrates that (N ₁ + N ₂)/P vs. I plot does not form a collapsed line and its pattern is similar to vertically flipped μ ₁ vs. I plot (Figure 4A). Therefore, (N ₁ + N ₂)/P may depend only on |σ _xz /σ _zz | or μ ₁.

If (N ₁ + N ₂)/P is indeed a function of |σ _xz /σ _zz | or μ ₁, we need one more equation to represent N ₁ and N ₂ separately. (N ₁ − N ₂)/P may provide that equation. However, Figures 9C, D show that (N ₁ − N ₂)/P plots whose horizontal axis is I or |σ _xz /σ _zz | do not have a well-collected collapse in contrast to the collapse seen in Figure 9B. (N ₁ − N ₂)/P is almost a constant between 0.1 and 0.15 except for the chute data which go up to 0.2. It seems as if the peculiarity of the chute data is canceled in the (N ₁ + N ₂)/P plot and pushed into the (N ₁ − N ₂)/P plot to become more conspicuous. Although (N ₁ − N ₂)/P cannot be written as a function of a single dimensionless variable, we can still utilize Θ _xx − Θ _yy in Figure 7A. For example, Figure 9E demonstrates that subtracting 1.6 ( Θ x x − Θ y y ) 1 / 6 from (N ₁ − N ₂)/P reduces the chute data peculiarity and collapses the data into a curve that depends only on I.

Using (N ₁ + N ₂)/P and (N ₁ − N ₂)/P data collapses, we may solve for μ ₂ and μ ₃ in terms of I, Θ and Θ _xx − Θ _yy . However, as mentioned above, Θ _xx − Θ _yy in a more complex flow geometry is not simple to define. Also, the data collapse in 9E is not yet clear because it is achieved only for two big branches of (N ₁ − N ₂)/P while adding an arbitrary function with two fitting parameters. Thus, further research is needed to build a more general rheological model that incorporates the temperature anisotropy and its impact on (N ₁ − N ₂)/P.