A mathematical model of the viscoelastic phenomenon employing derivatives of fractional order is examined in light of its consistency with thermodynamic principles. In particular, the development of constraints on parameters of the model ensure that the model predicts a nonnegative rate of energy dissipation and a nonnegative internal work. These constraints lead the model to predict realistic sinusoidal response as well as realistic relaxation and creep responses. Coupled with the established properties of the model, these attractive characteristics enhance its credibility for engineering analyses.