In this paper, rheological equations at different plastic conditions are derived on the basis of visco-elastoplastic Nishihara model, and then these equations are further implemented into discrete element method. By studying the numerical integral scheme of force-displacement relationship between particles, the methods for developing contact constitutive model in Particle Flow Code (PFC) are summarized. Thus a user defined visco-elastoplastic Nishihara contact constitutive model is proposed under PFC in 2 dimensions (PFC2D). Through stress relaxation experiments between two fixed balls, the accuracy of the developed model is verified for three cases. Numerical parameters are obtained by the uniaxial and triaxial compression experiments on rock specimens from underground roadway in Zhangjiawa Iron Mine, and then it is applied to simulate the uniaxial compressive creep experiment. The feasibility of Nishihara model in PFC2D for rock rheological properties is proved by comparing the results with experimental data, and then failure modes of a D-shaped roadway are investigated using Nishihara model and a linear-elastic model. Numerical results of Nishihara model are in good agreement with the theoretical result of Hoek-Brown strength criterion. The deformations of two sides are substantially greater than those of the top and bottom and the corresponding failures of two sides are more serious due to the creeping of wall rock. Numerical results, to large extent, agree with the field monitoring data. It is suggested that the support on both sides should be greatly strengthened for satisfying requirements of long-term stability.