  The dynamic damage model in viscoelastic materials in Kelvin-Voigt rheology
is discretised by a scheme which is coupled, suppresses spurious numerical
attenuation during vibrations, and has a variational structure with a convex
potential for small time-steps. In addition, this discretisation is numerically
stable and convergent for the time step going to zero. When combined with the
FEM spatial discretisation, 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 which are
computationally used in some engineering simulations in a non-reliable 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.
