Accurately simulating viscoelastic wave propagation in the subsurface is essential to image subsurface reservoirs. One established method to investigate constant-Q effects is based on the variable-order fractional viscoelastic wave equation. However, the mixed-domain operators introduced in this approach are cumbersome to calculate and impede a high-efficiency numerical solver. To solve this issue, we approximate the mixed-domain operator by using a polynomial approximation and derive the corresponding viscoelastic wave equation in anisotropic attenuative media. The constant-order fractional polynomial method is further developed to improve the approximation accuracy. Based on these polynomial approximations, only several constant-order Laplacians remain in the new form of the viscoelastic wave equations. Therefore, these equations can be easily solved using Fourier transform and finite-difference methods. In addition, we also perform an accuracy analysis among the conventional quadratic, cubic, and constant-order fractional polynomial approximations. It indicates that the relative error of the constant-order fractional polynomial method is much smaller than that generated by the conventional quadratic and cubic polynomial approximations. Compared to the variable-order fractional viscoelastic wave equation, our constant-order viscoelastic wave equation can simulate seismic wave propagations with comparable accuracy but with high efficiency. Numerical results based on our approach can well match reference solutions in homogeneous and heterogeneous models, demonstrating the accuracy and effectiveness of the new method.