A three-dimensional, incompressible, multiphase particle-in-cell method is presented for dense particle flows. The numerical technique solves the governing equations of the fluid phase using a continuum model and those of the particle phase using a Lagrangian model. Difficulties associated with calculating interparticle interactions for dense particle flows with volume fractions above 5% have been eliminated by mapping particle properties to an Eulerian grid and then mapping back computed stress tensors to particle positions. A subgrid particle, normal stress model for discrete particles which is robust and eliminates the need for an implicit calculation of the particle normal stress on the grid is presented. Interpolation operators and their properties are defined which provide compact support, are conservative, and provide fast solution for a large particle population. The solution scheme allows for distributions of types, sizes, and density of particles, with no numerical diffusion from the Lagrangian particle calculations. Particles are implicitly coupled to the fluid phase, and the fluid momentum and pressure equations are implicitly solved, which gives a robust solution.