Variational multiscale (VMS) methods were developed in the last decades and have been proven to work well for turbulent flows. In this work, we analyze the efficiency and accuracy of the residual-based VMS method and the orthogonal-subscales residual-based VMS method implemented in a matrix-free fashion with a sum-factorization approach. The methods are implemented in the open-source Computational Fluid Dynamics (CFD) software Lethe that solves the incompressible Navier-Stokes equations by applying a continuous Galerkin Finite Element discretization and it is built on top of the deal.II finite element library. Different benchmarks that consider wall-bounded flow, flow separation from a curved surface and energy dissipation are used for the validation. The results using different Reynolds numbers are compared to reference numerical and experimental data. The accuracy in terms of average quantities, such as velocities, Reynolds stresses and reattachment points is evaluated. In addition, a computational performance analysis is conducted in order to evaluate the scalability and parallel performance of the code when run in computer clusters to highlight the differences between the matrix-free and the traditional matrix-based approach. For this, a computational cost defined as wall-clock time times the number of cores utilized is reported in CPUh along with an efficiency that takes into account the discretization, the linear solver and the implementation. The metrics for the matrix-based approach are obtained using an implementation of the methods in the same code to allow for a fair comparison.