  In this paper, we design and analyze a Virtual Element discretization for the
steady motion of non-Newtonian, incompressible fluids. A specific
stabilization, tailored to mimic the monotonicity and boundedness properties of
the continuous operator, is introduced and theoretically investigated. The
proposed method has several appealing features, including the exact enforcement
of the divergence free condition and the possibility of making use of fully
general polygonal meshes. A complete well-posedness and convergence analysis of
the proposed method is presented under mild assumptions on the non-linear laws,
encompassing common examples such as the Carreau--Yasuda model. Numerical
experiments validating the theoretical bounds as well as demonstrating the
practical capabilities of the proposed formulation are presented.
