Using rigid body dynamics to solve multi-impact problems poses many difficulties and unanswered questions. Experiments and numerical simulation by using all kinds of compliant contact models clearly show that the outcomes of the post-velocities are not only influenced by the local dissipated energy, but also significantly effected by the coupling between various contacts. Based on the idea that the local dissipated energy depends on the constitutive relationship of the compliances, while the couplings among contacts exhibit the wave behaviors, this paper presents a new method that can well deal with multi-impact problems and produce energetically consistent and unique solutions to the post-impact velocities. Stronge's energetic coefficient is used as the energetic constraint to reflect the local dissipated energy at each impact, and the wave effects are coupled into rigid body models by using a distributing law that is associated with the relative contact stiffness and the relative potential energy stored at contact points. The evolution of energy is mapped into the impulsive-velocity level by stretching the time scale into the impulsive scale such that multi-impact processes can be described as a set of differential equations with respect to a normal impulse. Combining the distributing law with the Darboux-Keller's method of taking the normal impulse as an independent "time-like" variable can make multi-impact problems calculable and respect the energy constraints. The guidelines related to the selection for the independent normal impulse are presented and a numerical algorithm is developed in this paper. Comparisons between the theoretical predictions and the experimental results found in the literature as well as the numerical results obtained with compliant contact models are presented. Robustness analysis through numerical simulations is carried out and shows that the method proposed in this paper can well keep the essential mechanical features of the real system as it respects the evolution of motion in the energy sense. We use this method to investigate some interesting phenomena found in granular systems, and hope to extend it to the modeling of multiple impacts with friction, an issue that will be tackled in the future.