This paper considers viscoelastic incompressible solids in large deformations. A three-dimensional thermodynamical model is first proposed. It leads to a nonlinear variational problem which is approximated in space by mixed finite elements and in time by a first order implicit scheme. This approximation is analyzed on the corresponding linearized problem. Finally, an elimination strategy is described which reduces each step of the discrete problem to the solution of a standard hyperelastic problem to be solved by Newton's method.