An accurate and robust numerical method has been developed to simulate turbulent two-phase flows. The phase interface is tracked by the level-set method to capture frequent topological changes due to breaking or merging. Because of the broad-band characteristics of length scales in two-phase flow, a Lagrangian drop breakup model has been developed, which is coupled to the level-set method. In this approach, small subgrid droplets produced from resolved ligaments are then transferred from the level-set representation to the Lagrangian particles. The further secondary atomization is handled by a stochastic breakup model. When pinching-off of ligaments is not resolved on the level-set grid, a capillary breakup model is used to predict the drop size distribution from the pinching off and inserted as Lagrangian drops. This method improves the mass conservation as well as reducing the computational cost. For a high-fidelity simulation of two-phase flow, a new numerical algorithm has been developed to improve the robustness of the numerical method. The conservative formulation of Navier-Stokes equations is solved with a density correction term in the present method. The density flux terms are calculated from the level-set field for accuracy. In addition, a constant coefficient Poisson system is solved for pressure to satisfy the continuity equation in the fractional-step method. In order to show the capability of the method as an efficient tool in the breakup process, the atomization of a round liquid jet surrounded by a coaxial gas is considered. The numerical results are consistent with the observed breakup mechanisms in the experiment and the stability analysis. The drop size distribution of the resulting spray after breakup is also compared with the experimental data. The subgrid drops are also predicted by the Lagrangian drop breakup model, which shows the applicability of our method for numerical simulation of the atomization process. Both theoretical and numerical approaches are employed to investigate the stability mechanisms of the air layer drag reduction (ALDR) phenomenon. A linear viscous stability analysis is performed by solving the Orr-Sommerfeld equations in a two-dimensional two-phase Couette-Poiseuille flow configuration that mimics the far-downstream region from an air injector. Air-layer stability is reduced as the free-stream velocity, Froude number, and velocity gradients at the air-liquid interface are increased, whereas the air-layer stability is enhanced as the gas flow rate and surface tension force are increased. Nonlinear stability characteristics are also studied using numerical simulations with the same Couette flow configuration as indicated in the linear stability analysis. The study shows that the Weber number has a significant effect on the breakup of the phase interface. As the Weber number increases, the liquid ligaments become thinner, requiring higher grid resolution. Therefore, for simulations of high Weber number flows, the use of a Lagrangian spray breakup model is essential to predict the dynamics of subgrid-scale liquid structures. Direct Numerical Simulation (DNS) of two-phase flow is also performed to investigate the air layer drag reduction (ALDR) phenomenon in turbulent water flow over a backward-facing step. The Reynolds and Weber numbers based on the water properties and step height are 22,800 and 560, respectively. The total number of grid points is about 271 million for DNS. Two different air-flow injection rates are examined to investigate the mechanism and stability of the air layer. For high air-flow rate, the stable air layer is formed on the plate and more than 90% drag reduction is obtained, whereas, in the case of low air-flow rate, the air layer breaks up and ALDR is not achieved. The initial Kelvin-Helmholtz instability causes the streamwise wave structure, while turbulence interaction forms the spawise waves and causes ligament breakups. However, overall rupture of the air layer is mainly determined by the stability of the streamwise wave. The stability of the streamwise wave can be predicted from the stability analysis in the far-downstream region.