The aim of this work is to describe a numerical framework for reliably and robustly simulating the different kinematic conditions exhibited by granular materials while spreading ---from a stagnant condition, when the material is at rest, to a transition to granular flow, and back to a deposit profile. The gist of the employed modeling approach was already presented by the authors in a recent work (Cante et al., 2014), but no proper description of the underlying numerical techniques was provided therein. The present paper focuses precisely on the detailed discussion of such numerical techniques, as well as on its rigorous validation with the experimental results obtained by Lajeunesse, et al. in Ref. ( Lajeunesse et al., 2004). The constitutive model is based on the concepts of large strains plasticity. The yield surface is defined in terms of the Drucker Prager yield function, endowed with a deviatoric plastic flow and the elastic part by a hypoelastic model. The plastic flow condition is assumed nearly incompressible, so a u - p mixed formulation, with a stabilization of the pressure term via the Polynomial Pressure Projection (PPP), is employed. The numerical scheme takes as starting point the Particle Finite Element Method (PFEM) in which the spatial domain is continuously redefined by a different nodal reconnection, generated by a Delaunay triangulation. In contrast to classical PFEM approximations ( Idelsohn et al., 2004), in which the free boundary is obtained by a geometrical technique (a-shape method), in this work the boundary is treated as a material surface, and the boundary nodes are removed or inserted by means of an error function. One of the novelties of this work is the use of the so-called Impl-Ex hybrid integration technique to enhance the spectral properties of the algorithmic tangent moduli and thus reduce the number of iterations and robustness of the accompanying Newton-Raphson solution algorithm (compared with fully implicit schemes respectively). The new set of numerical tools implemented in the PFEM algorithm – including new discretization techniques, the use of a projection of the variables between meshes, and the constraint of the free-surface instead using classic a-shape – allows us to eliminate the negative Jacobians present during large deformation problems, which is one of the drawbacks in the simulation of granular flows. Finally, numerical results are compared with the experiments developed in Ref. (Lajeunesse et al., 2004), where a granular mass, initially confined in a cylindrical container, is suddenly allowed to spread by the sudden removal of the container. The study is carried out using different geometries with varying initial aspect ratios. The excellent agreement between computed and experimental results convincingly demonstrates the reliability of the model to reproduce different kinematic conditions in transient and stationary regimes. Peer Reviewed