The accurate and efficient modeling of granular flows and their interactions with external bodies is an open research problem. Continuum methods can be used to capture complexities neglected by terramechanics models without the computational expense of discrete element methods. Constitutive models and numerical solvers are the two primary aspects of the continuum methods. The viscoplastic size-dependent nonlocal granular fluidity (NGF) constitutive model has successfully provided a quantitative description of experimental flows in many different configurations in literature. This research develops a numerical approach, within a hyperelasticity framework, for implementing the dynamical form of NGF in three-dimensional material point method (3D MPM, an appropriate numerical solver for granular flow modeling). This approach is thermodynamically consistent to conserve energy, and the dynamical form includes the nonlocal effect of flow cessation. Excavation data, both quantitative measurements and qualitative visualization, are collected experimentally via our robotic equipment to evaluate the model with respect to the flow geometry as well as interaction forces. The results are further compared with the results from a recent modified plastic Drucker–Prager constitutive model, and in other configurations including wheel–soil interactions, a gravity-driven silo, and Taylor–Couette flow.