A three-dimensional granular contact dynamics scheme is presented. The scheme is variational in structure, thus making it possible to solve the governing equations by means of mathematical programming methods. To facilitate the modeling of natural grains using spherical geometries, a rolling resistance model is developed. A number of static and dynamic benchmark examples are considered including the granular column collapse problem where the agreement between simulation and previously published experimental results is found to be very good.