The challenge of predicting velocity and stress fields in any flowing granular material has proven to be a difficult one, from both computational and theoretical perspectives.  Indeed, researchers are still in search of the ``Navier-Stokes"-equivalent for flowing granular materials.   Granular flows can be adequately predicted using grain-by-grain discrete element methods (DEM), but these approaches become computationally unrealistic for large bodies of material and long times.   A robust continuum model, once identified, would have the practical benefit that it could be implemented at a meso-scale saving many orders of magnitude in computation time compared to DEM.

Here, we begin by synthesizing a 3D elasto-viscoplastic law for steady granular flow, merging an existing "frictional fluid" relation with a nonlinear granular elasticity relation to close the system.  The flow rate vanishes within a frictional (Drucker-Prager) yield surface and the elastic response is based on a mean-field model generalizing Hertz's contact law. The resulting form is general, able to produce flow and stress predictions  in any well-posed boundary value problem.  We implement it using ABAQUS/Explicit finite-element package and run test simulations in multiple geometries. The solutions are shown to compare favorably against a number of experiments and DEM simulations.

While this relation appears to function well for rapid flows, experimental results can often differ from the predictions in regions of slower flows.  We have been able to attribute many of these phenomena to nonlocal effects stemming from the finite-ness of the grain size.  To correct this, we consider the addition of a simple nonlocal term to the rheology in a fashion similar to recent nonlocal flow models in the emulsions community.  The results of this extended model are compared against many DEM steady-flow simulations in three different 2D geometries.  Quantitative agreement is found for all geometries and over various geometrical/loading parameters.  By natural extension, the nonlocal model is then converted to three dimensions with minimal changes, and is implemented numerically as a User-Element in the ABAQUS package.  We show that a single calibration of the 3D model quantitatively predicts hundreds of experimental flows in different geometries, including, for the first time, the wide-shear zones observed in the split-bottom annular couette cell, a geometry made infamous for resisting a theoretical treatment for almost a decade.