A finite-difference collocated scheme is proposed for the simulation of unsteady incompressible flows in cylindrical coordinates. The Rhie-Chow flux-interpolation technique is used to avoid the odd-even decoupling of the pressure field. It is shown that the calculation of the Coriolis and centrifugal acceleration terms based on the cell-center values lead to numerical instabilities near the pole in the inviscid regime. To avoid this issue, a new formulation based on the combination of cell-center and cell-face velocity components is proposed for the calculation of these terms. Furthermore, a new pole treatment methodology is introduced to calculate the viscous terms, resulting single-valued properties at the pole. As regards the extensions of the present scheme to higher orders of accuracy, a fourth-order spatial discretization is also introduced and verified. Moreover, it is shown that for higher-order extensions, the present pole treatment is able to preserve the order of accuracy and the L-2-norm and L-infinity-norm errors are greatly improved in comparison with the previous schemes. The robustness of the present method is assessed via five different numerical tests including laminar flows in a polar cavity, vortex-wall interaction, evolution of a monopole vortex-system to a tripole vortex-system, advection of the lamb dipole, and large eddy simulation (LES) of turbulent pipe flows. Compared to the previous LES results for turbulent pipe flows, significant improvement is obtained for the mean and fluctuating velocities especially at the core region of the channel. Similarly, at the wall region, the fluctuating axial and azimuthal velocities and the skewnesses of the axial and radial velocities are very close to the previous direct numerical simulation and experimental data. (C) 2015 Elsevier Inc. All rights reserved.