In this article, a novel constitutive model using fractional calculus is developed to capture the stress relaxation behavior of glassy polymers, where a variable‐order differential operator based on Marchaud fractional derivative is adopted. To assess the validity of the proposed model, a series of stress relaxation tests of the representative glassy polymer, poly(ethylene glycol‐co‐1,4‐cyclohexanedimethanol terephthalate), are conducted, and the stress responses of the specimens under different ambient temperatures are obtained. In view of the basic theory of fractional viscoelasticity, the varying order is assumed to bea function of time. Through the comparison of the fitting effect and reasonability of four kinds of order function based on the experimental data, it is concluded that the linear order function of time is suitable for our model which shows high accuracy. Moreover, the physical significance of order function is analytically derived, and the rising order is found to have a direct connection with the continuous softening of polymers during relaxation. Deeper investigations manifest that high ambient temperature accelerates the material softening, and the change of microstructures inside polymers could be reflected by the varying order, which is possible to provide a new perspective on manufacturing proper polymers with excellent resistance to stress relaxation.