We discuss a general purpose multibody solver which is able to simulate systems with more than 200,000 colliding rigid bodies with friction and joints. Multibody systems with rigid contacts entails the solution of nonsmooth equations of motion: the dynamics is non-smooth because of the
discontinuous nature of noninterpenetration, collision, and adhesion constraints. These lead to formulations based on non-linear complementarity problems: the more the contacts, the more difficult is the solution. Large number of contacts can occur, for instance, when simulating rock soil dynamics, conveyor belts, packaging devices, size-segregation shakers and granular flow problems which involve thousands or millions of particles in contact (ex: powder dynamics, tire-sand interaction, etc.).

An interesting case of granular flow is the refueling process in PBR nuclear reactors, where our method, based on rigid body dynamics, does not require small time steps as other methods (ex. Discrete Element Methods,DEM) which use spring-dashpot approximations. The huge amount of non-smooth
constraints arising from unilateral contacts with friction give rise to a NCP (Nonlinear Complementarity Problem) which we solve by means of a custom iterative method. Special care has been paid in implementing a high-performance matrix-less algorithm with O(n) space complexity. This solver has been implemented on multi-core systems and currently is under
porting to GP-GPU parallel architectures with stream processors. The proposed method has been implemented as a C++ API within the Chrono::Engine middleware.