Suspensions of rigid particles can be found in a wide range of natural and industrial settings. Detailed models of the movement of individual particles can be a useful tool to predict and understand macroscopic properties of these suspensions. Boundary integral equations have been shown to be an efficient method to do this, however for dense suspensions avoiding particle overlaps requires excessively small time steps. By applying a constraint requiring no contact at each time step to the variational Stokes equations we are able to maintain separation of all the particles without requiring a small time step. This enables long time horizon simulations of dense suspensions of rigid particles. After presenting the method and demonstrating its robustness we demonstrate its applicability to modeling physical phenomena such as effective viscosity and particle alignment angles.