Abstract: The Navier-Stokes equations are the fundamental equations in fluid mechanics. and Their numerical simulations are of very important practical significance. This paper proposes a three-step grad-div stabilization method for numerically solving the Navier-Stokes equations., which combines the advantages of both the grad-div stabilization method and the two-grid finite element discretization method, introducing the grad-div stabilization term to reduce the influence of pressure on the approximate velocity solution. Its main steps are as follows: firstly, by solving a grad-div stabilized nonlinear Navier-Stokes problem on a coarse mesh to get a rough solution, and then in following two steps, two grad-div stabilized and linearized Navier-Stokes problems are solved to obtain a finial solution on a fine grid. We estimate the error bounds of the approximate solution, and perform numerical experiments on known analytical solution examples and a forward step flow problem to verify the correctness of the theoretical analysis and the effectiveness of the proposed method.