The numerical method which is based on flux difference splitting, LU decomposition, and implicit high-resolution third-order Essentially Non-Oscillatory (ENO) scheme was constructed for the efficient computation of steady state solution to three-dimensional incompressible Navier-Stokes equations in general coordinates. The flowfields over underwater axisymmetric bodies, full-appended axisymmetric bodies and axisymetric bodies with a ring-wing duct were simulated. The method is proved to be capable of predicting the circumferential-mean velocity distribution at model scale to the accuracy of around 3% of measured values, and of predicting some details of flow features, for example, the wake harmonics.