We present a new family of algorithms for incompressible 3D Navier-Stokes equations in cylindrical geometry. A model problem of turbulent flow calculation in an infinite circular pipe {(r, ϕ, z): 0 ≤ r ≤ R, 0 ≤ ϕ < 2π, |z| < ∞} is considered and used for accuracy, stability, and efficiency estimations. Algorithms are based on Galerkin trigonometric approximation for uniform variables ϕ, z, on pseudospectral polynomial approximation in the r -direction (with different sets of collocation nodes) and on implicit and predictor-corrector time advancement schemes. In all cases high (infinite order) spatial accuracy is retained despite the presence of coordinate singularity at r = 0. To achieve this we exploit the behaviour of analytic functions of variables r, ϕ, z in the vicinity of r = 0. We analyze the advantages and disadvantages of four Navier-Stokes algorithms. In method A a new splitting technique is developed which makes use of a second-order predictor-corrector scheme and nontraditional fractional step procedure. Stability and efficiency characteristics of this scheme exceed that of the usually used mixed Adams-Bashforth/Crank-Nicolson time advancement. To minimize errors due to splitting, algorithm B is suggested that has no fractional steps. In this method pressure values are eliminated from discretized Navier-Stokes equations by means of equivalent matrix operations. Although conventional Chebyshev collocation nodes rl = R cos(πl/2Q), l = 0, 1, ..., Q, are used in both methods, the discrete boundary conditions at r = 0—consistent with analytic behaviour of solutions for small r—are fully accessible for the first time. In addition, approximations developed prevent the appearance of various pathological (with, e.g., spurious, parasitic modes, etc.) discretizations of Navier-Stokes operators. In algorithm C we propose a new set of collocation nodes rl = (1 - xl )R/2, l = 0, 1, ..., Q, where xl ϵ (-1, 1), l = 1, 2, ..., Q - 1, are the zeros of Jacobi polynomial P(2,1)0-1(x), xD = -1, x0 = 1. It is demonstrated that pseudospectral polynomial approximation with this set of nodes possesses the discrete analogue to the energy conservation law of the original Navier-Stokes initial boundary value problem. The latter is of special significance for the algorithm's nonlinear stability. In method D new dependent variables are introduced that completely consider the form of the analytic pressure and velocity components at small r. We show that the discrete Navier-Stokes equations admit in this case an efficient solution procedure. Finally, we present a technique that can be used for exhaustive a priori estimates of the algorithm's accuracy and stability characteristics in the linear approach.