In this work we present a variational derivation of a mixed formulation of the Kelvin- Voigt viscoelastic model and its numerical discretization. The Kelvin-Voigt is one of the most commonly used linear models of viscoelasticity and finds applications in many fields such as, e.g., propagation of shear waves in anisotropic and heterogeneous tissues such as kidney and skeletal muscles, dynamics of soil-structure interactions, or propagation of earthquake waves. Unlike the linear elasticity equation, the Kelvin-Voigt model cannot be casted as the Euler- Lagrange equations of an energy functional, and thus the standard duality technique, typically used to derive a mixed formulation, cannot be used directly. To circumvent this problem, we employ a technique recently developed for the solution of on an old De Giorgi’s conjecture, and recast our problem as the limit of a sequence of convex optimization problems. Then the aforementioned standard duality approach yields a sequence of saddle point problems whose limit is our final mixed formulation. The derived mixed formulation of viscoelasticity involves Hdiv regularity of the total stress (elastic and viscous) tensor together with H1 regularity of both the displacement and the velocity fields. This regularity is naturally called for by the physical principle of involving continuity of normal stresses across material interfaces. The numerical solution of this mixed dual formulation is obtained by low order finite element approximations for stresses, velocities and displacements. In order to avoid sat- isfying the inf-sup conditions required when the standard Galerkin approach is used, we employ a stabilized finite element formulation based on the variational multiscale concept that permits to use the same interpolation for all the unknowns. Time integration is performed using classical finite difference schemes.