In this research, using a phenomenological strain energy function, a 3D visco-hyperelastic constitutive equation for the finite deformation of elastomers is developed. For the quasi-static response, a strain energy function, composed of two exponential terms, is employed. The same energy function is considered as the kernel in the hereditary integral approach to derive the visco-hyperelastic model for the time-dependent response. To reduce the number of material parameters, an inverse power-law function between relaxation time and strain rate is proposed. Discretization of the analytical formulation of the visco-hyperelastic model in time is performed to be used in numerical simulations. This approach provides a recursive formulation to update the stress in each time step based on the deformation history. Then, the model is used to study the behavior of elastomeric bushing in radial, torsional and axial deformations. In all these cases, the material parameters are determined by applying a multi-objective optimization algorithm, in which the objective function is the difference between experimental and numerical results. The proposed visco-hyperelastic model shows accurate results for elastomers under 3D quasi-static and time-dependent loadings. Moreover, the accuracy of the results is not dependent on the amount of strain, the strain rate, and modes of deformation.