In this paper, we study the numerical solutions of viscoelastic bending wave equations$u_{t}(x,~t)-\int_{0}^{t}[\beta_{1}(t-s)\,u_{xx}(x,~s) - \beta_{2}(t-s)\,u_{xxxx}(x,~s)]ds = f(x,~t),$for $ 0<x<1,~ 0<t\leq T $, with self-adjoint boundary and initial value conditions, in which the functions $ \beta_{1}(t) $ and $ \beta_{2}(t) $ are completely monotonic on $ (0,~\infty) $ and locally integrable, but not constant. The equations are discretised in space by the finite difference method and in time by the Runge-Kutta convolution quadrature. The stability and convergence of the schemes are analyzed by the frequency domain and energy methods. Numerical experiments are provided to illustrate the accuracy and efficiency of the proposed schemes.