An abstract framework of the virtual element method is established for solving general elliptic hemivariational inequalities with or without constraint, and a unified a priori error analysis is given for both cases. Then, virtual element methods of arbitrary order are applied to solve three elliptic hemivariational inequalities arising in contact mechanics, and optimal order error estimates are shown with the linear virtual element solutions. Numerical simulation results are reported on several contact problems; in particular, the numerical convergence orders of the lowest order virtual element solutions are shown to be in good agreement with the theoretical predictions. This is a joint work with Bangmin Wu and Weimin Han.