A difference finite element method based on the mixed finite element pair \(((P_1^b,P_1^b,P_1) \times (P_1,P_1,P_1))\) - \((P_1 \times P_0)\) is presented for the three-dimensional stationary Navier–Stokes equations with damping. Moreover, based on this proposed method, a two-level discretization is constructed, which involves solving a problem of the Navier–Stokes equations with damping on coarse mesh with mesh sizes H and \(\mathcal {T}\) , and a general Stokes problem on fine mesh with mesh sizes \(h = O(H^2)\) and \(\tau = O(\mathcal {T}^2)\) . This two-level difference finite element method provides an approximate solution with the same convergence rate as the difference finite element solution, which involves solving a problem of the Navier–Stokes equations with damping on fine mesh with mesh sizes h and \(\tau \) . Hence, it can save a large amount of computational time. Finally, all computational results support the theoretical analysis and show the effectiveness of the two-level difference finite element method for solving the considered problem.
