A numerical modeling approach was developed to simulate the propagation of shear and longitudinal sound waves in arbitrary, dense dispersions of spherical particles. The scattering interactions were modeled with vector multipole functions and boundary condition solutions for each particle. Multiple scattering was simulated by translating the scattered wave fields from one particle to another with the use of translational addition theorems, summing the multiple-scattering contributions, and recalculating the scattering using an iterative method. The theory and initial results for the model are presented, including an integral derivation for the translational addition theorems. The model can simulate 3D material microstructures with a variety of particle size distributions, compositions, and volume fractions. To test the model, spectra and wave field patterns were generated from both ordered and disordered microstructures containing up to several hundred particles. The model predicts wave propagation phenomena such as refractive focusing, mode conversion, and band gap phenomena. The convergence of the iterations ranges from excellent to fair, and is dependent on the field (longitudinal or shear), particle configuration, and elastic wave frequency. The model is currently limited by the computation of sufficiently high multipole order for the simulation of dense particle dispersions.