This work is concerned with the numerical solution of the equations of the dynamics of flexible multibody systems with dissipative physical mechanisms. Specifically, we consider systems composed of rigid and deformable bodies made of a viscoelastic continuum, that may experience large displacements and strains, connected by joints. Within this framework, we present a novel integration scheme for these types of systems that intrinsically satisfies the laws of thermodynamics and existing symmetries. The resulting solutions are physically accurate since they preserve the fundamental physical properties of the model. Furthermore, the method gives excellent performance with respect to robustness and stability. Justification of these claims as well as numerical examples that illustrate the performance of the scheme are provided.