Abstract: The existing linear complementarity problem (LCP) modeling and numerical solution methods for non-smooth multibody system dynamics fail to accurately reflect the true physical characteristics. This study combines the Newton coefficient of restitution and the Moreau-Jean midpoint algorithm to develop the LCP based on the normal Lagrange’s equations. Meanwhile, the piecewise friction coefficient model is introduced to more realistically reflect the nonlinear characteristics of frictional forces. The results show that the LCP based on the normal Lagrange’s equations has superior numerical accuracy. The piecewise friction coefficient model enables this system to reach a stable state more quickly, in accordance with the energy dissipation law. This research is of great significance for understanding and predicting the dynamic behavior of multibody systems under complex contact conditions.