The flowing of powders brings a new challenging and interesting problem to the CFD community: at very high concentrations and low rate-of-strain, grains are in permanent contact, rolling on each other. Therefore a frictional stress model must be taken into account. This can be done using plasticity and similar theories in which the granular material behaviour is assumed to be independent of the velocity gradient or the rate-of-strain. This is in contrast to viscous Newtonian flow where the stress specifically depends on the rate-of-strain. Furthermore, unlike fluids, flowing powders do not exhibit viscosity and, again, this shows that a Newtonian rheology cannot accurately describe granular flow. It is assumed that the material is incompressible, dry, cohesionless, and perfectly rigid-plastic. Based on continuum theories,generalized Navier Stokes equations have been derived by replacing the velocity gradient by the shear rate, and the viscosity depends on pressure and shear rate. Thus the resulting equation is mathematically more complex than the Navier Stokes equations and is valid only when the material is deforming. In this note we present numerical algorithms to approximate these highly nonlinear equations. A Newton linearization technique is applied directly to the corresponding continuous variational formulations. The approximation of incompressible velocity fields is treated by using stabilized nonconforming Stokes elements and we use a Pressure Schur Complement smoother as defect correction inside of a direct multigrid approaches to solve the linear saddle-point problems with high numerical efficiency. The results of several computational experiments for realistic flow configurations are provided.