Computational fluid dynamics of viscous incompressible fluid flows is essential in many applied problems. Most often, we should focus on modeling turbulent regimes. Solutions at large Reynolds number values are of interest when a stationary solution of the boundary value problem for the Navier-Stokes equations exists, but this solution is unstable. In this paper, we propose iterative solution methods for such stationary solutions using natural pressure-velocity variables. These methods are not just theoretical but have practical applications. We consider two-dimensional problems for flows in a lid-driven cavity using finite element approximations on triangular irregular meshes. The developed computational algorithm is a powerful tool, allowing us to obtain two solutions when the Reynolds number exceeds the critical value for flows in a driven semi-elliptical cavity.