A comprehensive numerical model is developed using Lagrangian finite element (FE) formulation for investigating the steady-state viscoelastic (VE) rolling contact response. Schapery's nonlinear viscoelastic (NVE) model is adopted to simulate the VE behavior. The model accounts for large displacements and rotations. A spatially dependent incremental form of the VE constitutive equations is derived. The dependence on the history of the strain rate is expressed in terms of the spatial variation of the strain. The Lagrange multiplier approach is employed. The classical Coulomb's friction law is used. The developed model is verified and its applicability is demonstrated.