In this study, a new formulation of finite element method (FEM) has been extracted to analyze 2D viscoelastic problems. As there has not been enough accuracy and not sufficient literature in classical finite element modeling of viscoelastic problems, using a new set of shape functions founded on radial basis functions (RBFs) is recommended. Applying these new, RBF-based shape functions instead of the classical Lagrangian ones, results in subtler answers and conducts a reconsideration over the usual numerical method. Hankel functions are chosen, enriched and summed up with polynomial terms. Therefore, they satisfy not only polynomial terms, but also the first- and second-order Bessel functions simultaneously; which, in the case of classic shape functions, happens only for the polynomial function field. This method illustrates an approach with faster convergence rate and better robustness in different manners. Hence, it is less time-consuming and economical. Finally, various numerical examples are provided for the comparison of analytical solution, classic FEM and Hankel-based FEM, which show the much better agreement of the proposed method with analytical solution in comparison to classic FEM. Also, the number of nodes and degrees of freedom are reduced noticeably while maintaining accuracy in the interpolation of the adopted procedure.