With high order smooth non-uniform rational B-splines (NURBS) as basis function to simplify C¹ element construction, least squares isogeometric analysis is proposed for viscous incompressible Navier-Stokes equations. Governing equations are linearized by Picard or Newton method. Variational equation is derived from least squares functional defined by residuals of linearized equations. High order smooth finite dimensional spaces for velocity and pressure approximation are constructed by NURBS. Two benchmark flow problems were solved. Accurate numerical results were obtained for 2-dimensional lid driven flows. Global mass loss in flow past a cylinder in a channel decreased from 6% in classical least squares finite element method to 0.018%. It shows that the method is applicable to Navier-Stokes equations. It is better in mass conservation than least squares finite element method.