The principle objective of this dissertation was to develop a fast solution method for the accurate large scale modeling of fluid dynamics problems involving emulsions. The objective has been achieved by combining boundary element methods with fast iterative solvers based on the fast multipole method. Our solution method required extending the fast multipole method to generalized electrostatics problems to periodic domains. In addition, we developed an efficient hierarchical time stepping scheme and adopted accurate integration schemes, particularly useful in simulation of dense emulsions. Our solution method was implemented on parallel computers and used to model a creaming experiment of an emulsion.