A problem modelling fractional order viscoelasticity, is considered. Well-posedness of the model problem is studied within the framework of the semigroup theory of linear operators and the Galerkin approximation method. Spatial semidiscretization of the problem is formulated by means of finite element method, and full discretization of the problem is studied using continuous and discontinuous Galerkin methods. Stability estimates and optimal order a priori error estimates are proved.