The aim of this study was to develop an efficient and realistic numerical model in order to predict the dynamic response of belt drives. The belt was modeled as a planar beam element based on an absolute nodal coordinate formulation. A viscoelastic material was adopted for the belt and the corresponding damping and stiffness matrices were determined. The belt–pulley contact was formulated as a linear complementarity problem together with a penalty method. This made it possible for us to accurately predict the contact forces, including the stick and slip zones between the belt and the pulley. The belt-drive model was verified by comparing it with the available analytical solutions. A good agreement was found. Finally, the applicability of the method was demonstrated by considering non-steady belt-drive operating conditions.