Turbulent flows are governed by complex physical phenomena that have a wide range of scales in both time and space, making their simulation a challenging task even with modern high-performance computing architectures. In this work, we present a stabilized solver for the incompressible Navier-Stokes equations implemented in a matrix-free fashion with a sum-factorization approach. We solve the fully-coupled discretized problem in a monolithic way using a geometric multigrid preconditioner for the solution of the linear system. The solver is implemented in Lethe, an open-source Computational Fluid Dynamics (CFD) software, which uses the deal.II finite element library and a continuous Galerkin discretization. We verify the implementation using the method of manufactured solutions, demonstrating that the model recovers the appropriate order of convergence. In addition, we validate the results for two benchmarks, flow around a sphere and Taylor-Green vortex, to demonstrate the ability of the solver to accurately simulate problems that consider external flow and energy dissipation, at different Reynolds numbers, CFL and using different degrees for the equal-order velocity-pressure discretization. The results show better scalability and lower memory requirements than the ones of the matrix-based approach as the order increases.