A new approach, iteration projection method, is proposed to solve the saddle point problem obtained after the full discretization of the unsteady Navier-Stokes equations. The proposed method iterates projections in each time step with a proper convection form. We prove that the projection iterations converge with a certain parameter regime. Optimal iteration convergence can be achieved with the modulation of parameter values. This new method has several significant improvements over the Uzawa method and the projection method both theoretically and practically. First, when the iterative projections are fully convergent in each time step, the numerical velocity is weakly divergence free (pointwise divergence free in divergence free finite element methods), and the stability and error estimate are rigorously proven. With proper parameters, this method converges much faster than the Uzawa method. Second, numerical simulations show that with rather relaxed stopping criteria which require only a few iterations each time step, the numerical solution preserves stability and accuracy for high Reynolds numbers, where the convectional projection method would fail. Furthermore, this method retains the efficiency of the traditional projection method by decoupling the velocity and pressure fields, which splits the saddle point system into small elliptic problems. Three dimensional simulations with Taylor-Hood P2/P1 finite elements are presented to demonstrate the performance and efficiency of this method. More importantly, this method is a generic approach and thus there are many potential improvements and extensions of the iterative projection method, including utilization of various convection forms, association with stabilization techniques for high Reynolds numbers, applications to other saddle point problems.