The flow of a planar free shear layer with cylindrical bubbles is simulated using a finite difference/front tracking scheme. This approach allows direct numerical simulation of the multiphase flow by wholly incorporating the local bubble flow field in conjunction with the large scale vortical structures of the liquid. The role of large bubbles in modifying low Reynolds number (similar to 250) shear flow structures is investigated, specifically for bubbles whose diameter approaches the scale of the largest liquid eddies. The results indicate that duration of eddy crossing is the main mechanism for flow modulation, which is typically characterized by decreased vortex coherency and size, modified fluctuation statistics and significant variations in pairing/merging phenomena. The comparison of fluctuating statistics and flow field visualization also allowed qualitative discrimination between the modulation of the non-linear eddy dynamics and fluctuations due simply to the random bubble induced perturbations.