In this chapter, viscoelasticity effects in composites are studied. Three-dimensional CGs with linear viscoelastic matrices, containing linearly elastic spherical inclusions with or without interphases/coatings are treated. For each CG, the independent displacement fields are developed by the characteristic-length-scaled Papkovich-Neuber solutions and spherical harmonics. A compatible boundary displacement field is also assumed with Wachspress coordinates as nodal shape functions on each of the polygonal faces. Multi-field boundary variational principles are used to develop the CG stiffness matrices. After the establishment of CGs in Laplace transform domain, the homogenized and localized responses are transformed back to the time domain using the Zakian technique. 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 CGs are also used to study the effect of the negative Young’s modulus of particles on the stability and loss tangent of viscoelastic composites.