We consider a new Cartesian grid method for direct numerical simulations of fully coupled interaction of incompressible flow and spherical particles, based on a discontinuous extension of the pressure Poisson equation (PPE) across particle boundaries. We give a complete mathematical description of the boundary-integral treatment of the discontinuous PPE that includes the derivation of a new pressure boundary condition for accelerating boundaries and the solution of the system of boundary integral equations using spherical harmonics expansions. The model was validated with the standard test for finite Reynolds number flow around a sphere and with a novel test using the analytical solution for the Stokes flow past two adjacent spheres moving with the same velocity. The model capability and numerical efficiency was demonstrated with simulations for the collective settling of groups of 64–512 particles.