Abstract: A novel numerical algorithm of the elastoplastic damage constitutive model of rock based on the Mohr-Coulomb criterion is developed. The shortcomings of existing numerical algorithms are systematically analyzed in terms of whether the strain-softening process can be correctly simulated, whether the elastic modulus degradation caused by damage is considered, and the difficulty of numerical realization. The whole deformation and failure process of rock is divided into three stages, i.e., pre-peak plastic-strengthening stage, post-peak strain-softening stage and residual stage, and the corresponding numerical algorithm of each stage is deduced. For the pre-peak section, the fully implicit return mapping algorithm in the principal stress space is derived using three steps: elastic trial, plastic correction and damage correction. For the post-peak section, the plastic softening process is simplified into a series of brittle-plastic drop-plastic flow processes to overcome the problem of strain-softening simulation. The improved brittle-plastic method based on the plastic potential theory and the implicit algorithm for perfect plastic flow is introduced, and the elastic modulus degradation is considered by damage correction. For plastic flow in residual stage, it is simplified to two steps: elastic prediction and plastic correction. Through a secondary development of user subroutine UMAT of software ABAQUS based on FORTRAN language, the numerical simulation of elastoplastic damage constitutive model is realized. The correctness of the proposed numerical method is proven, as the triaxial compression cyclic loading and unloading test results and the analytical solutions of circular tunnel excavation considering strain-softening behavior agree well with simulation results. It is applied to the stability analysis of excavation and support of a large underground powerhouse groups to provide guidance for the design and construction.