We present a numerical method to deal efficiently with large numbers of particles in incompressible fluids. The interactions between particles and fluid are taken into account by a physically motivated ansatz based on locally defined drag forces. We demonstrate the validity of our approach by performing numerical simulations of sedimenting non-Brownian spheres in two spatial dimensions and compare our results with experiments. Our method reproduces qualitatively important aspects of the experimental findings, in particular the strong anisotropy of the hydrodynamic bulk self-diffusivities.