In this paper, we present a generalized model of a viscoelastic beam, which takes into account the influence of axial forces and incorporates a fractional constitutive relationship. In addition, we propose a novel numerical calculation method for analyzing fractional-order viscoelastic beams. This method takes the transverse displacement and bending moment of the beam as state variables, transforms the beam model into a discrete state space equation through the application of the central difference method, and utilizes an improved precise integration algorithm to solve the equation. To evaluate the performance of the method, the responses of the beam under two types of excitations, namely uniformly distributed transverse load and the motion of the support at both ends, are calculated under fixed hinge conditions. The results demonstrate that the present method has excellent accuracy and convergence, and also reveal some nonlinear phenomena of the system.