This paper aims to incorporate the fractional derivative viscoelastic model into a finite element analysis. Firstly, based on the constitutive equation of the fractional derivative three-parameter solid model (FTS), the constitutive equation is discretized by using the Grünwald–Letnikov definition of the fractional derivative, and the stress increment and strain increment relationship and Jacobian matrix are obtained by using the difference method. Subsequently, we degrade the model to establish stress increment and strain increment relationships and Jacobian matrices for the fractional derivative Kelvin model (FK) and fractional derivative Maxwell model (FM). Finally, we further degrade the fractional derivative viscoelastic model to derive stress increment and strain increment relationships and Jacobian matrices for a three-component solid model and Kelvin and Maxwell models. Based on these developments, a UMAT subroutine is implemented in ABAQUS 6.14 finite element software. Three different loading modes, including static load, dynamic load, and mobile load, are analyzed and calculated. The calculations primarily involve a convergence analysis, verification of numerical solutions, and comparative analysis of responses among different viscoelastic models.