Direct simulations of the interaction of a large number of deformable droplets are necessary for more accurate predictions of rheological properties and the microstructure of liquid-liquid systems. In the present study, a mathematical model of a three-dimensional flow of a mixture of two Newtonian liquids of a droplet structure in an unbounded domain at low Reynolds numbers is considered. An efficient computational method for simulation of the dynamics of a large number of deformable drops is developed and tested. This approach is based on the boundary element method for three-dimensional problems accelerated both via an advanced scalable algorithm (FMM), and via utilization of a heterogeneous computing architecture (multicore CPUs and graphics processors). This enables direct simulations of systems of tens of thousands of deformable droplets on PCs, which is confirmed by test and demo computations. The developed method can be used for solution of a wide range of problems related to the emulsion flow in micro- and nanoscales. Also it can be used to establish the closing relations for simulation of two-phase liquid-liquid flow based on the continuum approach in macroscales.