AbstractThis paper describes a method of performing the integration of generalized plasticity models, in which, unlike classical elastoplasticity, the yield surface is not explicitly defined. The algorithm is based on a generalized midpoint scheme and is applied to a specific generalized plasticity model for sands, in which a hyperelastic formulation is introduced to describe the reversible component of the soil response instead of the hypoelastic approach originally proposed. In this way, an efficient integration scheme is developed in the elastic strain space. The consistent, algorithmic tangent operator is derived. Isoerror maps are generated to study the local accuracy of the numerical integration algorithm. Results from a series of numerical examples based on the simulation of drained triaxial tests are given to illustrate the accuracy and convergence properties of the algorithm, both at the local and at the global level. Finally an example is given of the simulation of a cyclic triaxial test to illustrate the improvement on accuracy caused by the use of a hyperelastic law into the constitutive equations, as opposed to the hypoelastic formulation initially adopted in the model. Copyright © 2008 John Wiley & Sons, Ltd.