  In this paper, we introduce a new finite element method for solving the
Stokes equations in the primary velocity-pressure formulation. This method
employs $H(div)$ finite elements to approximate velocity, which leads to two
unique advantages: exact divergence free velocity field and
pressure-robustness. In addition, this method has a simple formulation without
any stabilizer or penalty term. Optimal-order error estimates are established
for the corresponding numerical approximation in various norms. Extensive
numerical investigations are conducted to test accuracy and robustness of the
method and to confirm the theory. The numerical examples cover low and high
order approximations up to the degree four, and 2D and 3D cases.

Numerical Analysis;Numerical Analysis