1 Introduction

Non-smooth dynamics is an important research field of multi-body system dynamics, and extensive works have been done for its modeling and numerical simulation from the past decade^[1-8]. Rigid and flexible multi-body systems with dry friction and impacts are non-smooth systems^[9]. Non-smooth events, like the state transitions of stick-slip between two bodies in the systems, make the dynamic equations of the systems discontinuous and nonlinear, complemented by inequality constraints due to unilateral contacts^[10]. A major difficulty in solving these equations results from instantaneous changes in the contact forces at transitions from slipping to sticking or transitions from sticking to slipping^[11]. For the wheeled multi-body systems, such as road vehicles, aircrafts and wheeled robots, when the wheels of the multi-body system roll on the ground, there are friction forces generated which are tangential to the surfaces of contact. The friction may lead to different states of motion, such as sticking (rolling without slipping, pure rolling) and slipping. When the wheels roll without slipping on the ground, the wheeled multi-body system is a multi-body system with bilateral nonholonomic constraints, and the wheels are subject to two-dimensional static friction. In order to study the dynamic control of these pure rolling wheeled multi-body systems, the wheeled multi-body systems were treated as multi-body systems with nonholonomic constraints, and extensive works had been done in recent papers^[12-15]. When the wheels roll with slipping, the system is a multi-body system with unilateral constraints, and the wheels are subject to two-dimensional kinetic friction. Tasora and Anitescu^[16] and Saux et al.^[17] studied the frictional contact problem of the wheeled multi-body systems and a rolling disk, respectively, where the rolling friction model was described as two-dimensional friction, and the rolling friction effect increased smoothly as the rolling speeds increased.

The friction with a stick-slip phenomenon is a major problem during the modeling and simulation of multi-body systems with friction, and it is important to choose a friction model capable of capturing its behavior and its influence on the system's dynamics^[18]. The Coulomb dry friction model is a fundamental and simple model of friction between dry contacting surfaces and can be used to describe the stick-slip phenomenon^[19]. Although the Coulomb dry friction model is straightforward, it presents some difficulties of implementation since it does not specify the friction force at the null velocity, and it is not easy to detect the stick-slip transitions in solving dynamic equations of the multi-body system. In order to avoid these difficulties, several friction models have been developed over the years based on modifications of Coulomb's approach^[20-22]. However, the modified Coulomb friction models cannot capture stick-slip motion. The Coulomb dry friction model will be used in the present paper.

Based on the Coulomb dry friction model, the problem of the transitions of the stick-slip of the wheels in the wheeled multi-body system is formulated as a linear complementarity problem (LCP) or a nonlinear complementarity problem (NCP) depending on the plane or spatial character of the contacts^[23]. When the friction is described as one-dimensional planar friction, the frictional contact problem can be solved by LCP methods^[24-25], but when the friction is described as two-dimensional friction, NCP methods are needed. However, solution methods for NCPs are still subjects of active research. Most methods use NCP-functions whose solution is zero, and others, like interior point methods, make use of penalty functions. In either case, solving NCPs is very cumbersome from the numerical point of view^[23]. In the present paper, by using the event-driven method, the problem of the transitions of the stick-slip of the wheels in the wheeled multi-body system is formulated as an NCP. The event-driven method can detect changes of stick-slip transitions, and resolve the exact transition time. Between the stick-slip events, the motion of the system is smooth and can be computed by a standard ODE-integrator with root-finding. If a stick-slip event occurs, the integration stops, and the computation of the contact forces is performed by solving an NCP^[26]. In order to solve the NCP, a linear approximation algorithm for the NCP is given, and the NCP will be solved by an iteration approach using the LCP method in the present paper.

This paper is organized as follows. In Section 2, in the wheeled multi-body systems, the wheel is treated as a rigid disk. The holonomic and nonholonomic constraint equations of the disk rolling without slipping on the ground are given as well as the correlations between Lagrangian multipliers and the constraint forces (holonomic constraint forces and nonholonomic constraint forces). Furthermore, the expressions of two-dimensional Coulomb dry friction are given when the disk is rolling with slipping. In Section 3, the dynamical equations of a multi-body system with one wheel are given. When the wheel slips, the dynamical equations of the multi-body system with holonomic constraints are obtained from the Lagrange equations of the first kind. When the wheel is in the pure rolling mode (rolling without slipping), the dynamical equations of the multi-body system with holonomic and nonholonomic constraints are obtained from the Routh equations. Some constraint-stabilized methods for multi-body dynamics with joints, contact and friction were studied in order to decrease the drift of the constraints^[27-28]. Baumgarte's stabilization method will be used to decrease the constraint drift. In Section 4, by using the event-driven method, the problem of the transitions of the stick-slip of the wheels in the wheeled multi-body system is formulated as an NCP. Then, a linear approximation algorithm for the NCP is given. The NCP will be solved by the iteration approach using the LCP method. In Section 5, a numerical example is considered as a demonstrative application example simulating a wheeled multi-body system consisting of a wheel on a shaft rotating around a pillar. The numerical results obtained show some dynamical behaviors of the wheeled multi-body system and constraint stabilization effects. Finally, the paper ends with some conclusions in Section 6.

2 Constraints and constraint forces

For a multi-body system with wheels, when a wheel rolls without slipping on the ground with dry friction, it is in the sticking state, and the static friction force acts on the wheel. When a wheel slips, it is in the slipping state, and the kinetic friction force acts on the wheel. In this paper, the wheels are treated as rigid disks, and the frictional contacts between wheels and the ground are characterized by a set-valued force law of the type of Coulomb's law for dry friction.

2.1 Holonomic and nonholonomic constraints

When a rigid disk rolls without slipping on the ground, the relative velocity of the disk point contact with the ground is zero. The disk is subject to holonomic and nonholonomic constraints.

Figure 1 shows that a rigid disk with the radius r rolls on the ground. An inertial coordinate frame Oxyz is attached to the ground. The position of disk mass center G is described by coordinates (x, y, z). As shown in Fig. 2, the body coordinate frame Gx'y'z' is fixed on the disk, which is obtained by rotating the frame GXYZ over angles (φ, θ, ψ). The orientation of the disk is determined using the following three steps:

Step 1  Rotate the frame GXYZ about the GZ axis through an angle φ, and the frame GP2Z is obtained.

Step 2  Rotate the frame GP2Z about the G2 axis through an angle θ, and the frame G123 is obtained.

Step 3  Rotate the frame G123 about the G3 axis through an angle ψ, and the frame Gx'y'z' is obtained.

When the disk purely rolls on the ground, the holonomic and nonholonomic constraints which the disk is subject to are as follows:

(1)

where f₁=0 is the holonomic constraint equation, which means that the disk remains contacting with the ground. The Lagrangian multiplier λ₁ is the constraint force corresponding to the equation f₁=0 and is the normal force acting on the contact point K of the disk shown in Fig. 1. f₁^* = 0 and f₂^* = 0 are the nonholonomic constraint equations to show that the disk rolls without slipping on the plane Oxy. The Lagrangian multipliers λ₁^* and λ₂^* are the constraint forces corresponding to equations f₁^* = 0 and f₂^* = 0, which are two components of the static friction force acting on the contact point K of the disk.

2.2 Two-dimensional Coulomb dry friction force and supporting force

When the rigid disk slips, the disk is subject to the kinetic friction force. According to the Coulomb dry friction law, the expressions of the friction force F_T and the supporting force F_N can be expressed as follows by introducing the coefficients of the static friction μ₀ and the kinetic friction μ:

(2)

where v_K is the velocity vector of the contact point K of the disk rolling on the plane Oxy, and .

When the disk rolls without slipping, it is in the sticking state. The friction force has a value within a range given by

(3)

v_K and can be expressed as and . v_K can be written as

(4)

When the friction force F_T can be expressed in the matrix form of

(5)

where F_T is the function of which can be expressed as the functions of (x, y, z, φ, θ, ψ) and

When the friction force F_T can be expressed in the matrix form of

(6)

where F_T is the function of , which can be expressed as the functions of When the problem of the stick-slip transitions of the contact point K can be converted into an NCP. Based on a linear approximation algorithm for the NCP, a numerical iteration method will be given to solve the NCP in Section 4 to determine the numerical results of and F_T.

When v_K = 0 and , the disk is subject to bilateral nonholonomic constraints. The friction force acting on the disk is the static friction force, which can be written as

(7)

3 Dynamic equations and algorithm

In the wheeled multi-body systems, the wheels have two motion states. One is the slipping state, and the other is the sticking state or pure rolling (rolling without slipping) state. When the wheel is in the slipping state, the dynamical equations of the multi-body system with holonomic constraints are obtained from the Lagrange equations of the first kind and a constraint stabilization method. When the wheel is in the pure rolling state, the dynamical equations of the multi-body system with holonomic and nonholonomic constraints are obtained from the Routh equations and a constraint stabilization method.

3.1 Dynamic equations and algorithm for slipping state

Let be the vector of generalized coordinates of the wheeled multi-body system with only one wheel. The wheel is treated as a rigid disk of the radius r. The configuration of the disk can be fully described by six independent variables , and other parts of the multi-body system are described by . Then, the vector of generalized coordinates can be written as .

Assume that the disk always keeps in contact with the ground. Then, the constraint equation of the disk can be written as

(8)

Let be the combined set of scleronomic holonomic constraints of the system besides the constraint equation of the disk f₁=0. Then, the scleronomic holonomic constraints of the system can be expressed together in the matrix notation,

(9)

Using the Lagrange equations of the first kind, the dynamic equations of multi-body system with holonomic constraints are expressed as

(10)

where is the column vector of generalized velocity of the multi-body system, is the Jacobian matrix of is the vector of Lagrange multipliers, in which represents the magnitude of the normal force on the contact point of the disk, and represents the constraint force of . F_T is the vector of the friction forces acting on the disk, the matrix W_T is a function of q, and the product of W_T and F_T represents the generalized forces of F_T. Q is the generalized force vector of active forces acting on the system.

L=T-V is the Lagrange function, in which T is the kinetic energy of the system, and V is the elastic potential energy of the system,

(11)

(12)

where is the generalized mass matrix, and is the generalized stiffness matrix.

Introducing Eqs. (11) and (12) into Eq. (10) leads to

(13)

where , and . To avoid the error accumulation due to constraint drift of the constraints, Baumgarte's stabilization technique is introduced,

(14)

where η and ξ are prescribed positive constants that represent the feedback control parameters for the velocity and position constraint violations^[29-31].

Equation (14) can be rewritten as

(15)

where .

Substituting Eq. (13) into Eq. (15) results in

(16)

By inserting λ₁=|F_N| into Eqs. (5) and (6), the two-dimensional Coulomb dry friction force F_T can be expressed as

(17)

(18)

From Eq. (4), v_K can be written as^[23]

(19)

where . Therefore, can be obtained from q and

From the above Eq. (19), can be written as^[23]

(20)

where . The computation of will be discussed in Section 4.

According to Eqs. (17) and (18), F_T can be written in the matrix form as

(21)

where the matrix A_FT is the function of q and .

Substituting Eq. (21) into Eq. (16) results in

(22)

For the sake of notation simplicity, Eq. (22) can be expressed in the matrix form as

(23)

(24)

(25)

Introducing Eq. (25) into Eq. (13) yields the differential equation of the multi-body system in the matrix notation,

(26)

Thus, the differential equation of motion can be solved by the numerical method for the ODE.

3.2 Dynamic equations and algorithm for sticking state

When the disk is in the pure rolling state , the multi-body system is subject to holonomic and nonholonomic constraints. With the Routh equations^[32],

(27)

where Φ^* is the column vector of nonholonomic constraints of the system. is the Jacobian matrix of is the vector of Lagrange multipliers about Φ^*, and it represents the constraint forces of Φ^*.

Equation (1) in Subsection 2.1 can be written as

(28)

Introducing Eqs. (11) and (12) into Eq. (27) leads to

(29)

Baumgarte's stabilization method can be used to decrease the constraint drift of the nonholonomic constraints, and this method can be written as

(30)

where γ is a prescribed positive constant that represents the feedback control parameter for the velocity constraint violations^[31], and .

Equation (30) can be rewritten as

(31)

where

Substituting Eq. (29) into Eqs. (15) and (31) respectively results in

(32)

(33)

Thus, Eqs. (32) and (33) can be written in the matrix form as

(34)

By introducing Eq. (34) into Eq. (29), the differential equations of the multi-body system with holonomic and nonholonomic constraints can be solved by the numerical method for ODEs.

4 Approximation algorithms for NCP

In the case of v_K = 0, there are two motion states when the wheel rolls on the ground, i.e., one is the slipping state and the other is the pure rolling state . The transitions of the motion states of the disk need to be detected first for solving the dynamic equations of the wheeled multi-body system numerically.

In this section, the two-dimensional Coulomb dry friction contact problem of multi-body systems with unilateral contacts is converted into an NCP.

4.1 NCP and two-dimensional Coulomb dry friction

The NCP is to find a vector such that

(35)

where F(x) is a mapping from R^k into R^k.

When , the vector , and the matrix , the NCP reduces to the LCP^[33].

The extended NCP^[34] is to find a vector such that

(36)

where G(x) and H(x) are mappings from R^k into R^k.

When v_K = 0, the conditions for the disk in the pure rolling state are that the friction force acting on the disk is less than the maximum static friction force, and can be written as .

The friction saturations are defined as^[26]

(37)

As shown in Fig. 3, F_T is the friction force vector, θ_K is the angle between F_T and the axis Ox, and the dashed line L coincides with F_T.

The components of F_T on the axes Ox and Oy can be expressed as

(38)

With projecting onto the line L, this projection value a_K can be solved by

(39)

a_K is splitted into the positive and negative parts^[35],

(40)

Thus, the friction saturations F⁺, F^- and a_K⁺, a_K^- are subject to the complementarity conditions,

(41)

According to Eqs. (13), (25), and (37)-(40), the friction saturations F⁺, F^- and a_K⁺, a_K^- can be expressed as

(42)

(43)

where H and G are vectors of 2 elements, and are general nonlinear functions of q and . Thus, the complementarity problem of two-dimensional friction belongs to the extended NCP as shown in Eq. (36). However, solving NCPs is very cumbersome from the numerical point of view^[23].

4.2 Linear approximation algorithm and iteration approach for NCP

When the disk rolls without slipping on the ground, v_K and are zero, and the magnitude of the static friction force can be described with the inequality which means that the values of the static friction force F_T are bounded by the circle with a radius of μ₀|F_N| as shown in Fig. 3.

The approximate approach to determining the value range of the static friction force is to replace the original circular zone of friction by an approximate polygonal zone^[36].

For example, a circular range can be replaced approximately by a square range as shown in Fig. 4 or an octagonal range as shown in Fig. 5. When the values of the static friction force F_T are bounded by the square with the side length of 2μ₀|F_N|, the value range of the static friction force can be expressed as

(44)

In Fig. 4, the x-and y-axes are two orthogonal axes. The inequalities in Eq. (44) mean that the magnitudes of the two orthogonal components of the static friction force are less than or equal to the magnitude of the maximum static friction force.

For each inequality in Eq. (44), the friction saturations are defined as^[32]

(45)

Equation (45) can be rewritten in the matrix notation,

(46)

By projecting onto the x-and y-axes, respectively, the values of the projection a₀ and a_π/2 are given by

(47)

a₀ and a_π/2 are splitted into the positive and negative parts^[33],

(48)

Thus, the vectors and are subject to the complementarity conditions,

(49)

In order to solve this complementarity problem, the correlations between should be clarified first.

Equation (48) can be written as

(50)

Thus, Eqs. (47) and (50) can be combined in the matrix form as

(51)

From Eqs. (46) and (51), the correlations between , as well as are obtained, respectively. However, the correlations between still need to be clarified.

Substituting Eq. (13) into Eq. (20) results in

(52)

However, Eq. (52) is still coupled with the variable λ.

For the sake of notation simplicity, Eq. (16) can be expressed in the matrix form as follows:

(53)

(54)

(55)

Substituting Eq. (55) into Eq. (52) yields

(56)

To reduce the above equation, it can be written as

(57)

where

Thus, from Eqs. (46), (51), and (57), the following functions can be obtained:

(58)

where A₁ and b₁ are the functions of q and

Thus, the square approximation algorithm is formulated as a four-dimensional linear complementarity problem (4D-LCP) by Eqs. (49) and (58).

Equation (58) can be solved by numerical methods such as Lemke's method^[24] to get the solutions and , and then , and can be obtained by solving Eqs. (50), (51), and (57), respectively.

If the magnitude of a₀ or a_π/2 is greater than zero, it means that the magnitude of is greater than zero, and the disk is in the slipping state.

If the magnitudes of a₀ and a_π/2 are both zero, it means that the values of the static friction force are bounded by the square with the side length of 2μ₀|F_N| as shown in Fig. 4. It is obvious that this approximation method has low accuracy.

When the values of the static friction force F_T are bounded by the octagonal range as shown in Fig. 5, this octagonal range is more approximate to the circular range than the square range as shown in Fig. 4.

Like the square approximation algorithm above, projecting F_T and onto the axes not only and but also the x-and y-axes as shown in Figs. 6 and 7, the vectors , and , and can be obtained by solving two complementarity equations of the 4D-LCP, respectively.

If all of the magnitudes of a₀, a_π/2, a_π/4, and a_3π/4 are zero, it means that the values of the static friction force are bounded by the octagon shown in Fig. 6. In this case, the disk can be treated as in the pure rolling state (sticking state).

When the values of the static friction force F_T are bounded by the regular polygon with 4n^* sides, and n^* is large enough, this regular polygon range is approximate to the circular range. Like the algorithm above, projecting F_T and onto the 2n^* axes which can be expressed by , j = 1, 2, 3, …, 2n^*, respectively (when , the axis is the y axis), can be obtained by solving n^* complementarity equations of the 4D-LCP. If all of the magnitudes of are zero, it means that the values of the static friction forces are bounded by the regular polygon with 4 × n^* sides, and the motion state of this disk can be treated as the pure rolling state (sticking state).

If some of the magnitudes of are greater than zero, it means that , and the disk is in the slipping state. In this case, there must be an with the largest magnitude among all the . When n^* is large enough, let and , then, can be solved by the following equation:

(59)

Using Eqs. (18) and (59), the friction force F_T acting on the disk can be computed.

4.3 Algorithm for multi-body system with unilateral non-smooth constraints and bilateral nonholonomic constraints

In Subsection 4.2, the stick-slip transition of a rolling disk on the frictional ground is formulated as LCPs.

The event-driven method will be used to detect changes of stick-slip transitions of the disk, and resolve the exact transition time. Between the stick-slip events, the motion of the system is smooth and can be computed by a standard ODE-integrator with root-finding. If a stick-slip event occurs, the integration stops, and the computation of the contact forces is performed by solving LCPs. Then, the dynamic equations of the wheeled multi-body system with two-dimensional Coulomb dry friction and nonholonomic constraints can be solved by the following steps:

Step 1  Give the geometric, inertial parameters of the system, and the initial conditions q₀, and initial n^*.

Step 2  Compute v_K to detect the motion state of the disk. If is sufficiently small), go to Step 3. If go to Step 4.

Step 3  When the disk is in the slipping state. Then, the friction force is given by Eq. (17), Eq. (26) is solved by the numerical method for the motion ODEs (like the Runge-Kutta method for an ODE) to determine q and at the next time, and go to Step 2.

Step 4  When , use the 4 × n^* polygonal algorithm in this paper to compute . If , go to Step 5. If , go to Step 6.

Step 5  When , check F_T computed by the 4 × n^* polygonal algorithm. If , then the disk is in the sticking state (rolling without slipping), Eq. (29) is solved by the numerical method for the ODE to determine q and at the next time, and go to Step 2. If , the values of the static friction force F_T are bounded by the regular polygon with 4 × n^* sides rather than the circle with the radius , in this case, the disk should slip, and should not be zero, then set n^* = n^* + 1, and go to Step 4 for computing and determining the friction force by Eq. (18).

Step 6  When , then the friction force is given by Eq. (18), Eq. (26) is solved by the numerical method for the ODE to determine q and at the next time, and go to Step 2. If , then set n^* = n^* + 1, and go to Step 4 for computing .

5 Numerical examples

The multi-body system with a rigid disk rolling on the ground is shown in Figs. 8 and 9. C₁ is the mass center of homogeneous disk. The disk of the radius r_C1 and the mass m_C1 is fixed on the massless round tube and rotates about the homogeneous shaft OA of the length l_A, the mass m_A and the diameter d_A. The disk C₁ and the shaft OA are connected by a spring of stiffness k and undeformed length l_S. The shaft OA hinges on a pillar Oz by the cardan joint. The coefficients of the kinetic friction and the static friction between the disk and the ground are μ and μ₀, respectively. The generalized coordinates of the system can be written as , where , and and are the mass center coordinates of the disk C₁ and shaft OA, respectively. and are the angles of the disk C₁ and shaft OA, respectively. The physical meaning of the generalized coordinates has been shown in Subsection 2.1.

The parameters of the system are given as follows:

The initial conditions of the system are given as follows:

The number of the polygon axes n^* = 8 in the polygonal approximation algorithm. When is the velocity magnitude of the disk contact point), is treated as zero. The feedback control parameters in the Baumgarte's stabilization technique η = 25, ξ = 50, and γ = 25.

The shaft OA is subject to a time-periodic torque M_OA which is parallel to the axis Oz and is expressed as

(60)

Case 1  The amplitude of time-periodic torque M_OA is M=5. The run time of the simulation is 10 s.

The motion of the disk is shown in Figs. 10 and 11 using the two methods. The solid line shows the motion of the disk using the approximation algorithm for an NCP. By contrast, the dashed line shows the motion of the disk using the trial and error algorithm. In the trial and error algorithm, the dynamic equations for the rolling disk are given using the Newton-Euler method.

Figure 10 gives the motion trajectory of disk mass center C₁, and Fig. 11 gives the time history of L_OC1 which is the distance between points O and C₁. The results obtained by the two methods are the same. When the value of L_OC1 is constant, the disk is in the pure rolling (sticking) state.

Figure 12 depicts the time history of the velocity magnitude of the disk contact point |v_K|. Figure 13 gives the time histories of the maximum static friction force μ₀|F_N| and the magnitude of friction force |F_T|. When the disk is in the sticking state. When , the stick-slip phenomena occur as shown in Figs. 12 and 13.

Case 2  The amplitude of time-periodic torque M_OA is M=25. The run time of the simulation is 25 s.

Figure 14 shows the motion trajectory of disk mass center C₁, and Fig. 15 illustrates the phase portrait of . Figure 16 shows the time history of L_OC1. Figure 17 gives the time histories of and . When the value of L_OC1 is constant, the disk is in the sticking state.

Figure 18 depicts the time history of the velocity magnitude of the disk contact point |v_K|. Figure 19 gives the time histories of the maximum static friction force μ₀|F_N| and the magnitude of friction force |F_T|. When , the disk rolls without slipping on the ground. When , the stick-slip phenomena occur.

These simulation results show that the disk has the periodic motion with stick-slip phenomena, and the periodic motion of the disk is the same as that of the excitation M_OA. During each period, the motion state of the disk changes three times from the pure rolling state (sticking state) to the slipping state.

Case 3  The amplitude of time-periodic torque M_OA is M= 3. The run time of the simulation is 50 s.

The disk is always in the pure rolling state when M=3. The drift of constraint equations using the two methods is shown in Fig. 20, where the dashed line shows the drift of constraint equations by the unstabilized method (η = 0, ξ = 0, γ = 0), and the solid line shows the drift by Baumgarte's stabilization method (η = 25, ξ = 50, γ = 25). Figure 20 illustrates that the drift of constraint equations can remain bounded with Baumgarte's stabilization method.

The motion trajectory of disk mass center is shown in Fig. 21 using two methods. The motion trajectory is a circle of the radius 1.5 m. Figure 22 shows the enlargement of the motion trajectory of disk mass center in the grey circle of Fig. 21. The dashed line shows that the radius of the disk trajectory changes using the unstabilized method, while the solid line shows that the radius of the disk trajectory remains constant using Baumgarte's stabilization method. Figures 21 and 22 illustrate the constraint stabilization effect of Baumgarte's stabilization method in this case.

6 Conclusions

In this paper, a numerical method for the dynamics of multi-body systems with two-dimensional Coulomb dry friction and nonholonomic constraints is presented based on the theory of the multi-body system dynamics and the algorithm for complementarity problems. The multi-body system with the wheels rolling on the ground with friction is considered. The frictional contacts between the wheel and the ground are characterized by a set-valued force law of the type of Coulomb's law for dry friction, and the friction belongs to two-dimensional friction. The constrains of the system include the constraints of perfect joints, the unilateral non-smooth constraints and the bilateral nonholonomic constraints between the wheel and the ground. The dynamical equations of the multi-body system are obtained using the Routh equations and a constraint stabilization method. Detecting the stick-slip transitions between the wheel and the ground is transformed into solving the NCP by the event-driven method. The complementarity conditions and complementarity equations of the two-dimensional friction problem are given. An iterative algorithm for solving the NCP is presented by the method for the LCP. Finally, a numerical example is considered as a demonstrative application example simulating a wheeled multi-body system consisting of a wheel on a shaft rotating around a pillar. The numerical results obtained show some dynamical behaviors of the wheeled multi-body system and constraint stabilization effects.