Sprays have wide applications in agriculture, pharmaceutical synthesis, engines, ink jet printing and so on. The successful spray applications and the control of spray param- eters require a thorough understanding towards the physical mechanisms. Numerical tools have been developed in the past few years for simulating the multiphase turbu- lent flows like sprays. Several researchers have successfully carried out direct numerical simulations (DNS) to investigate the primary breakup in such flows. DNS is accurate but requires extensive computational resources. In comparison, large eddy simulation (LES) is more practical, resolving only the large-scale flow structures and modeling the small-scale effects. The major difficulty with LES of multiphase turbulent flows is the need to model the interfacial subgrid-scale terms. Subgrid surface tension force, for ex- ample, plays an important role in the small droplet formation process. Subgrid surface tension force is, however, a highly non-linear term and can be difficult to model. In this research, we propose a new approach that combines the filtered density function (FDF) approach with the large eddy simulation. The major advantage of FDF is that the non-linear surface tension force appears in a closed form and thus needs no sub- grid modeling. The FDF transport equation is solved conveniently via a Lagrangian Monte-Carlo method. The Lagrangian approach is attractive in that it facilitates the transport of the liquid-gas interface without the diffusive or dispersive errors found in the Eulerian approaches. The surface tension source term in the momentum equation is closed using a Lagrangian volume of fluid (LVOF) approach. We utilize concepts from the smoothed particle hydrodynamics (SPH) in the LVOF approach to obtain the surface tension source term based on the Lagrangian particles. Several modifications have been made towards the original SPH formulation such that it is more suitable for the large-scale, turbulent multiphase flow simulations. Multiple particles are seeded in each Eulerian cell to achieve higher statistical accuracy, while the original SPH seeds one particle in each cell. What's more, a weighted SPH formula for the color function is adopted and is shown to be capable of handling variable particle number density. Performance assessment is via the rotation of Zalesak's disk and an oscillating elliptical droplet. Results show that the modified approach is able to handle the variable particle number density case appropriately. The simulations of multiphase turbulent flows are then performed with the proposed FDF-VOF methodology. At the same time, results from the simulations are compared with the DNS approach for validation and com- parison. Results show that the FDF-LES based approach can be a promising method, in that it models the flow with lower computational cost than DNS, yet maintaining accuracy in a model-free manor.