As a discontinuous deformation method, the block element method (BEM) characterizes a material’s elastoplastic behavior through the constitutive relation of thin-layer elements between adjacent blocks. To realistically simulate rock damage paths, this work improves the traditional BEM by using random Voronoi polygonal grids for discrete modeling. This approach mitigates the distortion of damage paths caused by regular grids through the randomness of the Voronoi grids. As the innovation of this work, the iterative algorithm is combined with polygonal geometric features so that the area–perimeter fractal dimension can be introduced to optimize random Voronoi grids. The iterative control index can effectively improve the geometric characteristics of the grid while maintaining the necessary randomness. On this basis, a constitutive relation model that considers both normal and tangential damage is proposed. The entire process from damage initiation to macroscopic fracture failure in rocks is described using two independent damage surfaces and a damage relationship based on geometric mapping relationships. The analysis results are in good agreement with existing experimental data. Furthermore, the sensitivity method is used to analyze the influence of key mechanical parameters in the constitutive model.