Summary: By combining grad-div stabilization used for improving pressure-robustness and full domain partition used for parallelization, a new parallel finite element method for the steady Navier-Stokes problem with friction boundary conditions is developed and analyzed. Within this parallel procedure, each processor assigned with one subdomain uses a composite grid to calculate a local approximation solution in its own subdomain, making the method simple and easy to carry out on the basis of existing sequential solver excluding a lot of effort in recoding on the top of existing serial software. We rigorously derive the uniform error estimates concerning the fine grid size \(h\), the viscosity \(\mu\) and stabilization parameter \(\alpha\) for the velocity, gradient of velocity and pressure in \(L^2\) norms and for the pressure in \(H^{-1}\) norm for the standard grad-div stabilized method. On the basis of the derived uniform error estimates and local a priori estimate for grad-div stabilized solution, we further give the uniform error estimates for the local and parallel grad-div stabilized methods. It is shown theoretically and numerically that our present method can reduce the influence of pressure on the error of velocity when the viscosity \(\mu\) is small. Specifically, it can yield better approximate velocity than its counterpart excluding grad-div stabilization. The smaller the viscosity coefficient, the higher improvment in accuracy of the approximate velocity. Besides, it can provide approximate solutions for the velocity and pressure with the same basically precision and convergence rate as the ones calculated by the standard grad-div stabilized method with a massive decrease in CPU time.