Abstract The dynamic damage model in viscoelastic materials in Kelvin–Voigt rheology is discretized by a scheme that is coupled, suppresses spurious numerical attenuation during vibrations and has a variational structure with a convex potential for small time steps. In addition, this discretization is numerically stable and convergent for the time step going to zero. When combined with a finite-element spatial discretization, it leads to an implementable scheme and to that iterative solvers (e.g., the Newton–Raphson) used for the nonlinear algebraic systems at each time level have guaranteed global convergence. Models that are computationally used in some engineering simulations in a nonreliable way are thus stabilized and theoretically justified in this viscoelastic rheology. In particular, this model and algorithm can be used in a reliable way for a dynamic fracture in the usual phase-field approximation.