Laboratory investigations of relaxation damage properties of high viscosity asphalt sand (HVAS) by uniaxial compression tests and modified generalized Maxwell model (GMM) to simulate viscoelastic characteristics coupling damage were carried out. A series of uniaxial compression relaxation tests were performed on HVAS specimens at different temperatures, loading rates, and constant levels of input strain. The results of the tests show that the peak point of relaxation modulus is highly influenced by the loading rate in the first half of an L-shaped curve, while the relaxation modulus is almost constant in the second half of the curve. It is suggested that for the HVAS relaxation tests, the temperature should be no less than −15°C. The GMM is used to determine the viscoelastic responses, the Weibull distribution function is used to characterize the damage of the HVAS and its evolution, and the modified GMM is a coupling of the two models. In this paper, the modified GMM is implemented through a secondary development with the USDFLD subroutine to analyze the relaxation damage process and improve the linear viscoelastic model in ABAQUS. Results show that the numerical method of coupling damage provides a better approximation of the test curve over almost the whole range. The results also show that the USDFLD subroutine can effectively predict the relaxation damage process of HVAS and can provide a theoretical support for crack control of asphalt pavements.