Nonlinear viscoelastic response of reinforced elastomers is modeled using a three-dimensional mixed
finite element method with a nonlocal pressure field. A general second-order unconditionally stable
exponential integrator based on a diagonal Padé approximation is developed and the Bergström–Boyce
nonlinear viscoelastic law is employed as a prototype model. An implicit finite element scheme with consistent
linearization is used and the novel integrator is successfully implemented. Finally, several viscoelastic
examples, including a study of the unit cell for a solid propellant, are solved to demonstrate the
computational algorithm and relevant underlying physics