In this talk, we focus on considering mathematical and numerical approaches to the motion of nonlinear viscoelastic springs with dynamic contact. This contact model is formulated by a nonlinear ordinary differential equation system and a pair of complementarity conditions. We propose three different numerical schemes in which each of them is combined by several numerical methods. The convergence of numerical trajectories is proved. All the three schemes are implemented to display numerical simulations.