Now I need to use the polar coordinate to simulate some problems. However, the centerline singularities caused some difficulties in my computations. So anyone could give me some help on this field.

P.s. I have refered some paper,such as [1] Numerical treatment of polar coordinate singlularities. J Comput. Phys. V.157,789-795 [2] Numercial treatment of cylindrical coordinate centerline singularities, IJCFD,20001,Vol15,251-263

Thanks a lot!

Since you use polar coordinate (r,\theta), the N-S equations should have asymptotic form in the limit r->0, derived using L'hospital rule. Usually, the center of a polar coordinate is located at symmetry of the flow, so you may use symmetric condition there.

As r->0 there will be no angular (theta) variation, i.e,

as r -> 0, d/d(theta) -> 0

The singularity will be removed when you apply this.

This can't be true. Suppose you had a shearing flow along the axis. Then the magnitude of the flow will depend on the direction from the axis (theta). Suppose you had a stagnation point at the axis. Then the direction of flow will vary with direction from the origin.

I must retract what I posted before. The resolution of this problem is not as simple as I wrote. Sometimes when the physical picture has symmetry then it can be used to remove the singularity, which by the way is purely artificial. Since the Jacobian becomes zero at r=0 the tranformation (to polar/spherical coordinates) is not valid there. In such a case is there any validity to the tranformed pde in polar/spherical coordinates at r=0 ?