A new method has been developed for the numerical simulation of finite size dispersed bubbly flows. Finite size bubbles are bigger than the numerical grid cell size, meaning that Lagrangian particle tracking (LPT) cannot be used, but they are yet too small to be accurately resolved by interface tracking (IT) methods. The need for such an approach arises mostly if large eddy simulation (LES) paradigm is used to resolve turbulence since it requires fine grids, and bubbles could easily cover several grid cells. The proposed method inherits features of IT and LPT, and features two-way coupling. The governing equations for the continuous phase are solved on Eulerian static grid. Bubbles are assumed to have a spherical shape, and a color function is used to monitor bubbles' position and to update the physical properties in every grid cell same way as one does in IT. However, the color function is not updated by solving an advection equation, but is rather re-initialized at every time step by the new bubbles' position calculated by solving Newton's second law of motion for each bubble. A new method to calculate the undisturbed liquid velocity required to model the hydrodynamic forces is also proposed in this paper. Two methods are proposed to impose the two-way coupling (i.e. bubbles feedback on the fluid). The first approach considers the bubbles as spherical rigid particles, and bubbles' velocities are imposed to the computational cells that lie inside the bubbles. In the second approach, a correction force is added inside each bubble to model the filtered scales due to usage of coarse grid resolution per bubble. The method is validated against terminal velocity of an air bubble rising in stagnant water to test the buoyancy and drag forces and its lateral motion in linear shear flow to check the lift force. The results are compared with simulations clone with IT Constrained Interpolation Profile (CIP) method to verify the accuracy of the two-way coupling. (C) 2014 Elsevier Ltd. All rights reserved