Predicting the mechanical response and damage evolution of elastomers under large deformation is of great significance in engineering applications. In this work, a finite element (FE) scheme is formulated and used to simulate rate-dependent damage in elastomers. While based on the theoretical model of Lavoie et al. (Extrem Mech Lett 8:114–124, 2016) and maintaining the key features such as kinetics of chain scission and polydispersity, the FE scheme presented here includes the consideration of finite compressibility. Both implicit and explicit algorithms are derived and implemented as user subroutines in ABAQUS. Validated against existing numerical results as well as experimental data on homogeneous deformation, the capability of the FE scheme to solve problems involving inhomogeneous deformation is further explored by simulating samples with pre-existing defects. The numerical results can successfully capture several interesting phenomena, such as crack blunting, stress reduction near defect caused by damage, and rate-dependent damage evolution. Good agreement is also found with experimental data on the strain field near a crack tip.