This paper applies a numerical method of polynomial function approximation to the numerical analysis of variable fractional order viscoelastic rotating beam. First, the governing equation of the viscoelastic rotating beam is established based on the variable fractional model of the viscoelastic material. Second, shifted Bernstein polynomials and Legendre polynomials are used as basis functions to approximate the governing equation and the original equation is converted to matrix product form. Based on the configuration method, the matrix equation is further transformed into algebraic equations and numerical solutions of the governing equation are obtained directly in the time domain. Finally, the efficiency of the proposed algorithm is proved by analyzing the numerical solutions of the displacement of rotating beam under different loads.