Discrete granular models are a natural choice when simulating dense suspensions, when the small distances between the particles lead to dominant contributions by lubrication and contact forces. In such case one can get rid of the costly resolution of Navier-Stokes equations, using closed form expressions for lubrication terms. However, those terms diverge when two hard spheres approach contact, and there are issues when integrating them directly with the finite precision of floating point calculations. In this paper, we introduce a visco-elasto-plastic interaction model for suspended spheres, which combines lubrication and elastic-frictional contact behaviour depending on surface roughness. An integration scheme is proposed for that model. Unlike earlier methods. the scheme enables an unconditionally stable time-integration of the interactions. The case of perfectly smooth spheres (null roughness), namely, is integrated correctly. The theoretical results are well reproduced in benchmark tests on two-sphere systems: one sphere sedimenting on one other and two spheres in a shear flow. From these benchmark tests, we propose phase diagrams showing the interplay between viscosity, roughness and stiffness. The second test case highlights the origin of non-reversibility particle trajectories. It is controlled by the particle roughness for rigid particles, and by the particle deformation when the capillary number is higher than the relative roughness. (C) 2020 Elsevier B.V. All rights reserved.