Brownian dynamics algorithms integrate Langevin equations numerically and allow to probe long time scales in simulations. A common requirement for such algorithms is that interactions in the system should vary little during an integration time step; therefore, computational efficiency worsens as the interactions become steeper. In the extreme case of hard-body interactions, standard numerical integrators become ill defined. Several approximate schemes have been invented to handle such cases, but little emphasis has been placed on testing the correctness of the integration scheme. Starting from the two-body Smoluchowski equation, the authors discuss a general method for the overdamped Brownian dynamics of hard spheres, recently developed by one of the authors. They test the accuracy of the algorithm and demonstrate its convergence for a number of analytically tractable test cases. © 2007 American Institute of Physics.