The solution of the Navier-Stokes equations (NSE) is a ubiquitous goal across a wide range of industrial and scientific fields. The main challenges arise in the convection-dominated regimes, where the meshsize should comply with the Kolmogorov scale to accurately capture the complex flow features. However, this requirement yields prohibitively high computational costs, and using coarser meshes is not a practical alternative, as it might result in inaccurate simulations with spurious numerical oscillations. To overcome this problem, turbulence model strategies were introduced to resolved large-scales eddies while modelling the unresolved subgrid-scales effects. The simplest yet physically consistent approach is the Smagorinsky model, which can be traced back to Joseph Smagorinsky and was conceived for the numerical simulation of turbulence in weather forecasting. Convergence analysis for the Smagorinsky model can be found in the literature for example in the context of finite element method, see in (Chacon Rebollo & Lewandowski 2014). In this work, we consider the Smagorinsky model for the Navier-Stokes equations within a Virtual Element framework, introduced in (Berone et al. 2025, arXiv.2510.03563) and based on the divergence-free VEM developed in (Beirao Da Veiga & Lovadina & Vacca 2017-2018). The numerical experiments for turbulent flows exploited the support of general meshes to preserve the isotropy of the Smagorinsky model while accurately capturing turbulent features through mesh refinement. Under the classical assumption of small data, we prove existence and uniqueness of a solution, and obtain the expected a priori error estimates for the Smagorinsky model.