Bulk viscosity effects on the collective flow harmonics in heavy ion collisions are investigated, on an event by event basis, using a newly developed 2+1 Lagrangian hydrodynamic code named v-USPhydro which implements the Smoothed Particle Hydrodynamics (SPH) algorithm for viscous hydrodynamics. A new formula for the bulk viscous corrections present in the distribution function at freeze-out is derived starting from the Boltzmann equation for multi-hadron species. Bulk viscosity is shown to enhance the collective flow Fourier coefficients from \(v_2(p_T)\) to \(v_5(p_T)\) when \(% p_{T}\sim 1-3\) GeV even when the bulk viscosity to entropy density ratio, \(% \zeta/s\), is significantly smaller than \(1/(4\pi)\).