Abstract: Presently, dynamics of nonsmooth multibody systems is a hot research topic.The usual approach in treating such systems is to derive basic system fromthe original system by removing the nonsmooth constraints firstly. TheLagrange equations of the second kind of basic system combine with thecomplementarity condition of the nonsmooth constraints to set up at eachdiscrete moment in time a Linear Complementarity Problem (LCP).This article focuses on the problem of dynamic modeling and numericalsimulating of multibody systems with friction contacts. By neglecting theclearance and the effect of impact between rigid bodies and constraints, thestate variables in the differential equations are continuous. Due to theset-value mapping characteristic of dry friction forces, the differentialequations of motion have discontinuous right-hand vector fields, thereforeto allow our system to be classified as a Filippov system. In addition tothe friction constraints, our model also incorporates frictionless bilateralones. If the simulation has only friction constraints, then the problem isan ODE-LCP model. Combining friction contacts with frictionless bilateralconstaints, the ODE-LCP model has to be extended to DAE-LCP (DAE,Differential Algebra Equation) mixed model. In order to obtain the model,the basic system is derived from the original system by removing thefriction constraints firstly. Because the equality constraints are retainedin the basic system, the dynamic model of basic system is a set of DAE. Withthe aid of constraint Jacobian matrix, the normal contact forces andtangential friction forces of nonsmooth constraints, which obey thecomplementarity contact laws, are added to the DAE of the basic system toobtain the DAE-LCP mixed model.Approaches used in the past for simulating rigid multibody dynamics withfriction contacts include piecewise DAE approaches, acceleration-force LCPapproaches, and velocity-impulse LCP-based time-stepping methods.Recognizing that the nature of the frictional constraint can inducesick-slip motion, the last approach is used in this work, which has theadvantage that it does not suffer from the detection for stick-sliptransition that could appear in the first two approaches. This framework isbased on a LCP, but it is different from acceleration-force LCP approachesthat attempt to find the accelerations of the bodies.Our approach considersimpulses and velocities as the fundamental unknowns. Acceleration-force LCPapproaches solve for accelerations from the dynamics equations and then usethe accelerations in an integration procedure. Because the complementaritylaw between acceleration and friction surplus is valid only when therelative tangential velocity is zero, zero crossing detection for velocityis required. However the complementarity law between velocity and frictionsurplus remains valid through a full-range of motion, so in contrast toacceleration-force schemes, the velocity-impulse methods need noevent-detection. In the new framework, the integration and dynamicalresolution steps are combined. The main achievement of this approach is thatit has solutions for any configuration. As the time-step tends to zero, asubsequence of the numerical solutions approaches the solution of adifferential inclusion.Our method is carried out in a numerical example, and the simulation resultsindicate that this method is effective.