SUMMARYAn Eulerian-Lagrangian approach is developed for the simulation of turbulent bubbly flows in complex systems. The liquid phase is treated as a continuum and the Navier-Stokes equations are solved in an unstructured grid, finite volume framework for turbulent flows. The dynamics of the disperse phase is modeled in a Lagrangian frame and includes models for the motion of each individual bubble, bubble size variations due to the local pressure changes, and interactions among the bubbles and with boundaries. The bubble growth/collapse is modeled by the Rayleigh-Plesset (RP) equation. Three modeling approaches are considered: (a) one-way coupling, where the influence of the bubble on the fluid flow is neglected, (b) two-way coupling, where the momentum-exchange between the fluid and the bubbles is modeled, and (c) volumetric coupling, where the volumetric displacement of the fluid by the bubble motion and the momentum-exchange are modeled. A novel adaptive time-stepping scheme based on stability-analysis of the non-linear bubble dynamics equations is developed. The numerical approach is verified for various single bubble test cases to show second-order accuracy. Interactions of multiple bubbles with vortical flows are simulated to study the effectiveness of the volumetric coupling approach in predicting the flow features observed experimentally. Finally, the numerical approach is used to perform a large-eddy simulation in two configurations: (i) flow over a cavity to predict small-scale cavitation and inception and (ii) a rising dense bubble plume in a stationary water column. The results show good predictive capability of the numerical algorithm in capturing complex flow features.