Contact of viscoelastic materials with complicated properties and surface topography require numerical solution approaches. This paper presents a 3-D semianalytical contact model for viscoelastic materials. With the hereditary integral operator and elastic-viscoelastic correspondence principle, surface displacement is expressed in terms of viscoelastic creep compliance and contact pressure distribution history in the course of a contact process. Through discretizing the contact equations in both spatial and temporal dimensions, a numerical algorithm based on the robust Conjugate Gradient method and Fast Fourier transform has been developed to solve the normal approach, contact pressure, and real contact area simultaneously. The transient contact analysis in the time domain is computationally expensive. The fast Fourier transform algorithm can help reduce the computation cost significantly. The comparisons of the new numerical results with an analytical viscoelastic contact solution for Maxwell materials and with an indentation test measurement reported in the literature has validated and demonstrated the accuracy of the proposed model. Moreover, the present model has been used to simulate the contact between a polymethyl methacrylate (PMMA) substrate and a rigid sphere driven by step, ramped, and harmonic normal loads. The validated model and numerical method can successfully compute the viscoelastic contact responses of polymer-based materials with time-dependent properties and surface roughness subjected to complicated loading profiles.