We focus on the continuum simulation of dense granular flows, of which we note two characteristics: first, due to the Coulombic nature of the friction between particles, the material exhibits a yield-stress which depends of the local pressure. Then, contrary to ordinary fluids, granular flows can expend; however despite this observation, most previous approaches enforce a divergence-free velocity field. This has severe drawbacks, as it can lead to negative pressure values in the wake of obstacles, and therefore difficulties in the definition of the yield-stress. In this work we allow for positive values of the divergence. Numerically, most attempts at simulating granular media relied either on regularization techniques, or on an Augmented Lagrangian algorithm. We introduce a novel method based on nonsmooth optimization theory. Thanks to the absence of any regularization, we can accurately capture dead-zones and the transition between regimes.