Abstract: In the numerical solution of steady incompressible Navier-Stokes equations with frictional boundary conditions, the accuracy of velocity approximations obtained by the standard mixed finite element method depends on pressure. When the pressure is high and the viscosity coefficient is small, the accuracy of velocity approximations becomes unsatisfactory. To address this issue, this paper proposes a gradient–divergence (grad–div) stabilized two-level finite element method. The method employs a hierarchical strategy: first, a nonlinear Navier-Stokes-type variational inequality problem with a grad–div stabilization term is solved on a coarse grid to obtain a coarse-grid solution; then, the nonlinear convective term is linearized using the coarse-grid solution on a fine grid, and a linearized problem (Ossen, Stokes, or Newton type) is solved to derive the final solution. The grad–div stabilization term is utilized to suppress the influence of pressure on velocity. Theoretical error bounds for the approximate solutions of this numerical method are derived, and numerical experiments validate the correctness of theoretical results and the effectiveness of the method. Numerical results show that compared with the standard two-level finite element method, the proposed method significantly improves the accuracy of velocity approximations.