Particulate composites with polymer matrix are used in engineering applications, such as solid propellants and road pavements. The durability of such materials depend critically on the characterization of mechanical property degradation. Mechanical damage study of these composites thus become important. A phenomenological model describing the mechanical damage in polymer matrix is developed in the present work. Polymer is modelled as a hyperelastic - viscoelastic(VE) solid. The damage is characterized by the progressive degradation of storage modulus. This hyperelastic - VE model, with two term Prony series approximation of the relaxation modulus, is coupled with isotropic and anisotropic damage formulations. Two approaches of anisotropic damage evolution, viz. explicit and implicit, are implemented accounting thermodynamic considerations. Implicit isotropic and anisotropic damage evolutions are derived by a dissipation potential function based on the effective elastic stored energy density. An exponential evolution obtained from experiments is adapted in explicit approach and anisotropy induced by Poisson type coupling is made evident. The analytical VE - damage framework is numerically implemented as stress integration algorithms, and the corresponding three dimensional results are illustrated. The effect of including multiple terms in Prony series approximation and cyclic loading are also presented. The developed numerical model is employed to solve boundary value problems in commercial finite element analysis software ABAQUS® through user material subroutines (VUMAT). It is seen from the results that, the damage implicitly affects the relaxation time of the polymer. The stress hardening behaviour occuring immediately after the damage initiation is captured and subsequent softening is also demonstrated.