A physically based, isostructural, constitutive model is described for simulating the onset of nonlinear viscoelasticity in multiaxial creep of glassy polymers, as needed in stress analyses of load-bearing components. In the linear viscoelastic limit, shear response reduces to that of a generalized Maxwell model, while hydrostatic response is Hookean. Nonlinearity enters through Eyring-type rate process kinetics. The equations of the model are solved numerically using a pseudo-linear approximation through each time step, leading to an incremental equation for stress that would be convenient for use in finite element analyses. The model and its assumptions were tested using tension, shear and combined tension/shear creep experiments on well-aged poly(methyl methacrylate) at 70 degrees C. Reproducibility tests confirmed the assumption of constant glass structure for strains up to similar to 1.5 x 10(-2). Shear and pressure activation volumes were obtained by fitting the dependence of the shear compliance on stress invariants. The data showed unequivocally that shear activation volumes vary with log(relaxation time), and excellent agreement was obtained for a linear variation, consistent with the "compensation rule" of polymer thermo-viscoelasticity. The activation volumes are large (many monomer units), indicating the cooperative nature of the elementary flow process. Interestingly, they are of the same order as those applying to yield and plastic flow. Although the model finds success in simulating creep, it fails to describe so accurately the strain recovery on unloading Possible explanations are suggested.