Abstract: An implicit version of discrete element method (DEM) called non-smooth contact dynamics (NSCD) is used to simulate the crushing process of single rockfill particle under one-dimensional (1D) compression. The angular and irregular shapes of rockfill particles are represented using polyhedrons, which are discretized into smaller elementary cells through the Voronoi tessellation. The interaction between neighboring elementary cells is described by the cohesive zone model (CZM), where hybrid tensile and shear failure modes are considered. Consequently, particle crushing can be captured by the breakage of CZM bonds subjected to external loading. The proposed method for modeling the particle crushing avoids drawbacks of the traditional fragment replacement method and bonded sphere method in DEM. The Brazilian splitting tests on granite are performed to calibrate the parameters for CZM. 1D compression tests on single particles with different sizes and shapes further reveal that the crushing strength of particles follows the Weibull distribution, and both the magnitudes and the variances decrease with the increasing particle size. For the particles with different principal axes (e.g., elongated and platy), the crushing strength loaded in the major axis is the smallest among all directions. It is also shown that the elongated and platy particles have smaller average strength than the spherical ones given the same equivalent particle size.