The variable-order fractional (VOF) operator has been extensively employed to depict complex physical phenomena of viscoelastic materials, during which the reasonable form of the proper variable-order (VO) function will significantly affect the applicability and accuracy of proposed fractional model. To this regard, this work aims to provide a technical guidance for VO function selection in VOF viscoelastic modeling. In this paper, we develop a VOF model based on the classical Zener model, the suitable VO function form and the applicability of VOF model are investigated by means of characterizing the stress relaxation and creep recovery behaviors, in which the monotonicity of VO functions (decreasing or increasing) are discussed in detail.The results demonstrate that the stress relaxation and creep recovery responses have a realistic physical meaning only when the VOF rheological model with a monotonically decreasing VO function. Furthermore, three definitions of VOFD are also considered to determine the optimal VO function form and the reasonability of the proposed model is validated by fitting the experimental data of polymer matric materials. The analysis found that the decreasing exponential VO function exhibits good prospects in describing the time-dependent behaviors, and it was found that the V2 and V3 definition types performed better. The results also show that the proposed model has unique advantage compared to the other existing viscoelastic models.Finally, the proposed VOF model is employed to predict the viscoelastic behavior to test the model's predictive ability.