In this article, the 3D viscoelastic computational grains (CGs) with spherical inclusions, with and without interphases/coatings are developed to study the viscoelastic behavior of polymer-based heterogeneous materials.Papkovich-Neuber solutions and spherical harmonics are adopted to develop the independent displacement fields inside the grains, and Wachspress coordinate is used to interpolate compatible displacement fields at inter-element surfaces. The elemental stiffness matrices are developed with multi-field boundary variational principles in the Laplace domain, after which Zakian technique is adopted to invert both the homogenized and localized responses back to the time domain with good accuracy and efficiency. With different kinds of models to describe the property of the viscoelastic polymers, the generated homogenized moduli and localized stress distributions are validated against the experimental data, simulations by commercial FE software, and predictions by composite spherical assemblage models. Parametric studies are also carried out to investigate the influence of material and geometric parameters on the behavior of viscoelastic composites. Finally, the viscoelastic CG is also used to study the effect of the negative Young's modulus of particles on the stability and loss tangent of viscoelastic composites.