  Cylindrical coordinates are often used in computational fluid dynamics, in
particular, when one considers gas flow accreting onto a central object.
Although the cylindrical coordinates have several advantages in describing
rotation, they have apparent singularity along the axis at the coordinate
origin (z-axis). This singularity introduces difficulties in numerical
simulations. First, it is difficult to reproduce the flow across the z-axis.
Second, the time step is extremely shortened by the CFL condition near the
z-axis because the numerical cell thereof is narrow in the azimuthal direction
for a given angular resolution. Here, we propose a new discretization scheme to
overcome these difficulties. In our new scheme, we consider changes in the
direction of the unit vector within a cell when evaluating the flux across each
cell surface. Besides, we evaluate the source term in the radial component of
the momentum equation from the thermal and dynamic pressures working on the
azimuthal cell surface. The new scheme is designed to be free-stream preserving
so that flow with uniform density, pressure, and velocity isan exact solution
of the discretized equation. These improvements are essential to using a lower
angular resolution in the innermost area and thus to elongate each time step.
Our examples demonstrate that the innermost circular region around the axis can
be resolved by only six numerical cells. We present an application to an
accreting compact star surrounded by a disk, in addition to Sod shock tube and
rotating outflow tests.
