Using the discrete element method (DEM) with clustering, a novel means of numerically modeling damage of particles is presented. Damage, such as grain crushing, is treated by allowing clusters to break apart according to a failure criterion based upon sliding work. If the accumulated work done on an individual DEM particle of a cluster exceeds a threshold, that particle is allowed to break from the cluster. A value for the critical energy density is determined by comparing the degree of particle breakage from numerical simulations to data from laboratory tests. Numerical simulations were also conducted to determine the impact of particle damage on interface behavior. It was found that a very distinct shear zone was evident when particle damage was considered and that this occurred without significant reduction of the maximum shear strength of the medium. Also, the degree of damage was shown to be related to the angularity of the clusters.