Discrete element method (DEM), such as particle flow code (PFC), has become one of the current highlights to investigate the particle breakage process in the geotechnical field. In this study, based on the point load crushing criterion considering the local stress concentration, a three-dimensional (3D) particle crushing model was established. Particually, the particle balance among particles before and after crushing is ensured by the Apollo filling and expansion method, and a size effect factor is used to characterize strength of particles with different sizes. Then, numerical simulations were performed on three types of sand with different degrees of breakage, namely silica sand, calcareous sand and Sacramento River sand. The comparisons of numerical and experimental results show that the established model can well describe compressive properties of granular materials with different degrees of breakage. Comparing with the failure criterion rooted in Mohr–Coulomb theory on the basis of mean stress, the criterion in the established model shows a better manner in reflecting the real fragmentation phenomena under the same loading condition. At the same time, the established model reveals the influence of crushing on the anisotropic dissipation and gradation evolution of granular materials.