With the development of geothermal energy and unconventional reservoirs, shales have become the popular target formation. This kind of organic-rich or clay-rich geomaterial has significant viscoelastic behavior which may lead to severe accidents, due to the creep rupture and the relaxation of pre-stress. In order to understand the viscoelastic behavior and simulate it quickly and accurately, we have developed a novel procedure for solving viscoelastic problems based on fractional calculus. Fractional viscoelastic constitutive equations were built by introducing a springpot into the classical component models. Fractional and classical constitutive equations were compared by the fitting results of triaxial experimental data of shale samples. Fitting efficiency of fractional constitutive equation is enhanced if compared to classical equations. Besides, we also found that fractional order is directly correlated to the amount of total volume content of clay and kerogen. This implies that the physical meaning of fractional model is clearer and its extrapolation is more accurate. Fracture width creep model for unpropped fractures and closure stress relaxation model for well propped fractures were developed and analyzed. Their simulation results in terms of long-term and short-term scales differ tremendously and hence should be studied individually. The influence of lateral compliances of real shale samples on fracture width creep or closure stress relaxation is not significant and may be ignored in long-term simulations. Simulation results of closure stress relaxation could explain the abnormal proppant production, which occurs several days later than the start of natural flow-back.