Please tell me your idea and considerations in doing calculation on 3D Spherical coordinate by finite difference method

Everything is similar EXCEPT at the origin r=0, the difference method at r=0 is different and needs special treatment. This is the first thing I thought. Try to solve a simple heat conduction equation you may find the problem.

HL

When I make it into conservative form, In momentem equations, the velocity component are left behind. Can I treat them as production term?

Please check the paper "An accurate second-order approximation factorization method for time-dependent incompressible Navier-Stokes equations in spherical polar coordinates, W.sha, K.Nakabayashi and H.Ueda, J. Comput. phys.142, 47-66, 1998". Hope it would help you.

The N-S equations in the cylinderial and spherical systems can not be written into fully conservation form. I think you can treat the extra terms as production terms (source term).