The Non-Smooth Problem of Wheel Motion

Wiesław Grzesikiewicz¹, Artur Zbiciak^2*

¹Institute of Vehicles, Faculty of Automotive and Construction Machinery Engineering, Warsaw University of Technology, Warsaw, Poland.
²Institute of Roads and Bridges, Faculty of Civil Engineering, Warsaw University of Technology, Warsaw, Poland.

DOI: 10.4236/apm.2020.1011041                

Abstract

Planar motion of a non-deformable wheel under the action of non-ideal unilateral constraints is considered. The mathematical description of this phenomenon has a form of a non-smooth initial value problem. The non-smoothness of this problem means that its solution is determined by an absolutely continuous function having a discontinuous first derivative. For this reason, a collision problem describing abrupt changes of velocity has been formulated next to the equations of motion specifying the acceleration. The non-idealness of constraints means that the constraint reaction force includes also a component resulting from the friction between the wheel and the constraints. Differential equations specifying acceleration of the wheel making contact with the constraints and algebraic equations for determining the changes in the wheel’s velocity at the moment of collision have been formulated in the paper. The principal task in these formulations is to determine the reaction forces of the considered constraints. This task is specified by the relationships between acceleration and the constraint reaction force components. In the description of the collision, these relations refer to the post-collision velocities and reaction force impulses. For determining an approximate solution of the formulated wheel motion problem, an original numerical method and a computer program for wheel motion simulation have been developed. Selected results illustrating the changes in displacements and velocity have been presented.

Keywords

Non-Smooth Mechanics, Non-Ideal Constraints, Impact, Friction, Rigid-Body Dynamics

1. Introduction

In the classical (Newtonian) problem of the motion of a mechanical system with N degrees of freedom, we find—within a fixed time interval—the solution to the differential equation,

$\ddot{X}=f\left(t,X,\stackrel{\dot{}}{X}\right),t\in \left[0,t_{end}\right)$, (1a)

satisfying the initial conditions,

$X\left(0\right)=X_{\text{o}},\stackrel{\dot{}}{X}\left(0\right)=V_{\text{o}}$. (1a)

In the initial problem so formulated, it is assumed that the function $f:R^3\to R^N$ is continuous and that it satisfies the Lipschitz condition for the second and third arguments. Then, there exists a conclusive solution describing the motion of the system and the function X is continuous and double differentiable, that is, $X\in C^2\left(\left[0,t_{end}\right),R^N\right)$.

Non-smoothness of the problem of motion implies that the function f is discontinuous or non-differentiable with respect to its second or third argument. A problem of this kind applies, for instance, to mechanical systems with unilateral constraints or with dry friction forces. In the monograph [1] where a system with unilateral constraints is analysed it is demonstrated that the solution to the problem exists within the class of absolutely continuous functions $X\in C_{ab}\left(\left[t_{\text{o}},t_{end}\right),R^N\right)$. This means that the derivative $\stackrel{\dot{}}{X}$ is discontinuous, and, an additional task specifying discontinuous changes in velocity needs to be formulated.

In the further considered non-smooth wheel motion problem, we investigate non-ideal unilateral constraints. The said non-idealness implies that the reaction force of constraints is combined with the dry friction force perpendicular to the reaction force. A hypothesis of this kind was put forward at the close of the 19th century by Painlevé [2] [3]. The contemporary analysis of this problem is based on the fundamentals of the non-smooth analysis presented in works by Jean, Moreau and Panagiotopoulos [4] [5] [6] [7] [8], and in works by the authors of this paper on machine and vehicle modelling, for example [9] [10] [11]. One should also mention contemporary studies [12] [13] [14] in which a singular problem known as the Painlevé paradox, widely discussed at the beginning of the 20^th century by Prandtl, Hamel, Klein and Mises [2] is considered.

In relation to the above-mentioned studies, in this study we consider a specific method of description of motion with non-ideal unilateral constraints, suited for analysing a relatively simple problem relating to planar motion of a rigid wheel. The problem defined in this way fully illustrates the non-smooth issues of unilateral constraints while allowing leaving out a situation in which the Painlevé paradox occurs.

In addition to formulating the motion problem, we propose a method for determining the solution to this problem and present the results of wheel motion simulations.

2. Geometric and Kinematic Characteristics of the System

We consider a non-deformable wheel whose motion in the gravitational field is constrained by non-deformable planes with rough surfaces. Figure 1 is an illustration of the analysed system, with three coordinates describing the motion of the wheel, that is: $\left(x,y\right)$ —the coordinates of the wheel’s centre point, S and $\phi$ —the rotation angle of the wheel.

The motion constraints in Figure 1 are shown by the bold line representing the boundary of the area of the admissible positions of the wheel on the plane $O{x}_{\text{o}}{y}_{\text{o}}$. Figure 1 shows also the main dimensions of the boundary line and of the wheel, i.e. $r,h,l$, given that we are considering the case, when $r>h$ and $l>r$. The gravitational field vector g is also marked.

It follows from the preliminary analysis of Figure 1 that the constraints under consideration can be specified by the area of the admissible positions of the centre point of the wheel, since the third coordinate, i.e. the wheel rotation angle $\phi$, is not constrained.

Based on the geometrical analysis of Figure 1 we determine the description of the said constraint in the form of the relation

$\left(x\left(t\right),y\left(t\right)\right)\in \Omega \subset R^2\forall t\in \left[0,t_{end}\right)$, (2)

where the set $\Omega$, specifying the area of the admissible positions of the centre point of the wheel has the form

$\Omega :=\left\{\left(x,y\right):x\ge r;y\ge \Phi \left(x\right)\right\}$, (3)

if the function $\Phi$ is specified as

$\Phi \left(x\right):=\{\begin{array}{ll}h+r\hfill & \text{for}r\le x\le l_1\hfill \\ f\left(x\right)\hfill & \text{for}{l}_1<x\le l_1+a\hfill \\ r\hfill & \text{for}x>l_1+a\hfill \end{array}$, (4a)

where

$l_1:=l-r$, $a=\sqrt{r^2-{\left(r-h\right)}^2}$ (4b)

$f\left(x\right):=\sqrt{r^2-{\left(x-l_1\right)}^2}+h$. (4c)

Figure 2 shows the set of points $\Omega$ with the interior of the set $\Omega$ marked in grey, i.e. $\text{Int}\Omega$, and its boundary marked with the dashed line $\text{Fr}\Omega$. The

Figure 1. Geometrical representation of the considered wheel–constraints system.

descriptions of these sets have been established based on Equations (3) and (4), i.e.

$\text{Int}\Omega :=\left\{\left(x,y\right):x>r;y>\Phi \left(x\right)\right\}$ (5a)

$\text{Fr}\Omega :=\left\{\left(x,y\right):\left\{y\ge h+r,x=r\right\}\cup \left\{y=\Phi \left(x\right),x>r\right\}\right\}$. (5b)

Furthermore, Figure 2 shows also the positions of the wheel in selected sectors of the boundary $\text{Fr}\Omega$. A more detailed graph of wheel positions in all the characteristic sectors of the boundary $\text{Fr}\Omega$ are shown in Figure 3.

For the description of a constraint to motion to be complete, apart from the set $\Omega$, we need to establish also kinematic implications, which follow from the relation (2) and will be referred to as its consequences. In the study [3] the author considered in detail the consequences, which determine the sets of admissible velocities and admissible accelerations.

If the wheel (cf. Figure 2) is in such a position that $\left(x,y\right)\in \text{Int}\Omega$, i.e. when the wheel is not in contact with the boundary of a constraint, then the above-mentioned kinematic constraints do not occur, and the wheel is moving freely.

Figure 2. Graph of the set $\Omega$ and its boundary $\text{Fr}\Omega$ and a number of specific positions in which the wheel is in contact with the boundary of constraint.

Figure 3. Six specific positions in which the wheel is in contact with the boundary of constraint.

When the wheel makes contact with the boundary of a constraint, i.e. $\left(x,y\right)\in \text{Fr}\Omega$, the velocity of the centre point of the wheel $\left(v_x,v_y\right)$ is constrained, as it should satisfy the relation

$\left(v_x,v_y\right)\in D\Omega \left(x,y\right)\subset R^2,\left(x,y\right)\in \text{Fr}\Omega$, (6)

where $D\Omega$ is the mapping referred to as the first consequence of relation (2). The definition of the mapping is given in [1], while a detailed description of this mapping in the problem under consideration will be provided further in this paper.

If the position and velocity of the centre point of the wheel satisfy the relation (6), then the acceleration of the centre point of the wheel can also be constrained. However, if

$\left(x,y\right)\in \text{Fr}\Omega ,\left(v_x,v_y\right)\in \text{Int}D\Omega \left(x,y\right)$, (7)

then such a constraint does not arise, as the wheel making contact with the boundary of constraints moves away from it.

A limitation to the acceleration of the centre point of the wheel arises when the wheel makes contact with the boundary of the constraint and the velocity of the centre point of the wheel belongs to the boundary of the set of constraints to velocity $D\Omega$, i.e. when

$\left(x,y\right)\in \text{Fr}\Omega \text{and}\left(v_x,v_y\right)\in \text{Fr}D\Omega \left(x,y\right)$. (8a)

In such a situation the wheel making contact with the boundary of constraint can either roll on it or remain still. Then the acceleration of the centre point of the wheel $\left(a_x,a_y\right)$ is constrained, as there should be satisfied the relation

$\left(a_x,a_y\right)\in D^2\Omega \left(\left(x,y\right),\left(v_x,v_y\right)\right)$, (8b)

where $\left(a_x,a_y\right)$ is the mapping referred to as the second consequence of the relation (2), determining the set of the admissible accelerations of the centre point of the wheel. A definition of this mapping is provided in [1] and a detailed description referring to the wheel will be given further in Section 3.

In the situation specified by Equation (8a), the constraints are acting, which means that the reaction force, which will describe further on, is arising out of them. During the motion of the wheel there can also occur a situation when the position and the velocity of the wheel are such that

$\left(x,y\right)\in \text{Fr}\Omega \text{,}\left(v_x,v_y\right)\notin D\Omega \left(x,y\right)$, (9)

which means that the velocity of the wheel is not admissible due to constraints, as it does not satisfy the condition (6). In this state, the wheel collides with rigid constraints. During the collision the velocity of the wheel changes sharply (abruptly) and reaches a value, which satisfies the relation (6). The formula describing the change in velocity during the collision will be presented further on in the paper.

As mentioned before, the consequences $D\Omega$ and $D^2\Omega$ specify the sets of admissible velocities and accelerations of the wheel making contact with the constraining surface (Figure 1). The principles of determining these mappings have been described in detail in [1] [9]. Taking into account the relatively simple form of the boundary of the set $\Omega$ (Figure 2), we will determine the descriptions of these mappings in the form which requires no explanation.

3. Model of the Action of Constraints on the Wheel

An action of a constraint on the wheel occurs when the wheel comes into contact with the constraint. In this position, there can be generated a reaction force acting on the wheel, which we will discuss in detail further on in the paper.

Figure 4 is an example illustration of a wheel making contact with a plane constraining its motion at point A. A local coordinate system Ant, whose axis An is perpendicular, and axis At is tangent to the boundary $\text{Fr}\Omega$, is marked in the figure; also the velocity of the point A situated on the circumference of the wheel is shown in it. Due to the constraints to the motion of the wheel, the normal component of this velocity should satisfy the condition

$v_n\ge 0$. (10)

The velocity for which $v_n<0$ is in admissible due to constraints; in this situation a collision occurs, which will be considered further on in this paper.

If $v_n>0$, then the constraints do not limit the motion of the wheel, as the wheel moves away from the boundary.

Figure 4(b) shows the vector of acceleration of the point A.

If the velocity of point A is such, that $v_n=0$, then the constraint forces are acting on this acceleration

$a_n\ge 0\text{for}{v}_n=0$. (11)

In the above-described situation, there arises the reaction force illustrated in Figure 4(c). In the case of ideal constraints only contact force $R_n\ge 0$ is acting, while the tangential component, $R_t=0$ is not generated.

For non-ideal constraints a hypothesis is assumed according to which the tangential component maps the force of friction between the wheel and the boundary of constraints. The non-ideal constraints model so specified is described by relations, which determine the relationship between the reaction force R and the acceleration and velocity of the point A belonging to the wheel:

Figure 4. Velocity, acceleration and reaction vectors acting on the wheel making contact with the boundary of the set $\Omega$ at the point A.

- for the normal component $R_n$ we have:

$R_n\ge 0,a_n\ge 0,a_n{R}_n=0\text{for}{v}_n=0$, (12a)

- while the tangential component $R_t$, when $R_n>0$, is specified by the relations:

$R_t=\mu R_n\tau ,\tau :=\text{sign}{v}_t\text{for}{v}_t\ne 0$, (12b)

$R_t=\mu R_n\tau ,\tau \in \{\begin{array}{l}\left\{\text{sign}{a}_t\right\}\text{for}{a}_t\ne 0\\ \left[-1,+1\right]\text{for}{a}_t=0\end{array}\text{for}{v}_t=0$, (12c)

where $\mu$ —the coefficient of friction between the wheel and the boundary of constraints; $\tau$ —the friction multiplier.

The above relations can be illustrated with graphs; Figure 5(a) represents a graph of the relation (12a), while the relation (12c) is shown in Figure 5(b). Furthermore, the relation (12a) (Figure 5(a)) can be equivalently described by the relationship:

$R_n={\left[R_n-\rho a_n\right]}^{+}\text{,}\rho >0$, (13a)

where ${[\cdot ]}^{+}$ denotes such a function, that

${\left[\xi \right]}^{+}:=\{\begin{array}{ll}\xi \hfill & \text{for}\xi \ge 0\hfill \\ 0\hfill & \text{for}\xi <0\hfill \end{array}$. (13b)

Likewise, for the relation (11c) (Figure 5(b)) we obtain

$\tau =\Pi \left(\tau +\rho a_t\right)\text{,}\rho >0$, (14a)

where the function $\Pi$ has been defined as:

$\Pi \left(\xi \right):=\{\begin{array}{ll}\xi \hfill & \text{for}\left|\xi \right|\le 1\hfill \\ \text{sign}\xi \hfill & \text{for}\left|\xi \right|>1\hfill \end{array}$. (14b)

In the above-presented hypothesis on the interaction between non-deformable bodies in contact we have not taken into consideration the moment of forces, referred to as rolling resistances of the wheel. This is the moment of friction, the limit value of which is proportional to the contact force $R_n$. The direction of this moment is opposite to the angular velocity of the wheel. In order to simplify the description of the reaction force acting on the wheel, and assuming that the value of rolling resistance is very small, we have disregarded their effect on the motion of the wheel.

Figure 5. Graphs of relations (12a) and (12c).

Now we will consider a situation when the wheel collides with the constraints, that is when $v_n<0$. Schematic descriptions of the wheel’s velocities immediately before the collision $V^{-}$ and after the collision $V^{+}$ are shown in Figure 6(a) and Figure 6(b), while Figure 6(c) represents the reaction force impulses, ${\tilde{R}}_n$ and ${\tilde{R}}_t$, inducing the change in the velocity of the wheel. The friction impulse ${\tilde{R}}_t$ is determined according to the Routh’s hypothesis [3].

The already mentioned effect of the collision is described by a relation between the post-collision velocities $V^{+}$ (Figure 6(b)) and the reaction impulses (Figure 6(c)). This relation is formulated in accordance with the assumed hypothesis. For instance, when smooth (ideal) constraints are analyzed, the tangential component, that is ${\tilde{R}}_t=0$ is disregarded.

In this paper, we will present a description of the hypothesis on plastic collision, and on elastic collision or elasto-plastic collision, based on the assumption that the constraints are non-ideal.

A description of the plastic collision has a form analogical to the description of the reaction force, as given in the formulae (12). The plastic collision is specified by the relations:

${\tilde{R}}_n\ge 0,v_n^{+}\ge 0,v_n^{+}{\tilde{R}}_n=0$, (15a)

${\tilde{R}}_t=\mu {\tilde{R}}_n\tau ,\tau \in \{\begin{array}{ll}\left\{\text{sign}{v}_t^{+}\right\}\hfill & \text{for}{v}_t^{+}\ne 0\hfill \\ \left[-1,+1\right]\hfill & \text{for}{v}_t^{+}=0\hfill \end{array}$. (15b)

The relation (15b) has been determined based on the Routh’s hypothesis [3].

Graphs illustrating these relations are analogical to those shown in Figure 5 and they can be presented in an equivalent form, analogically to the formulae (13) and (14).

The description of the plastic or elastic-plastic collision is based on the Newton’s-Poisson’s hypothesis analyzed in [1]. According to this hypothesis the resultant reaction force impulse during such a collision is specified by the formulae:

${\tilde{R}}_n^{\ast }=\left(1+\beta \right){\tilde{R}}_n$, (16a)

${\tilde{R}}_t^{\ast }=\mu {\tilde{R}}_n^{\ast }{\tau }^{\ast }$, (16b)

${\tau }^{\ast }\in \{\begin{array}{ll}\left\{\text{sign}{v}_t^{+}\right\}\hfill & \text{for}{v}_t^{+}\ne 0\hfill \\ \left[-1,+1\right]\hfill & \text{for}{v}_t^{+}=0\hfill \end{array}$, (16c)

Figure 6. Velocities and impulses acting on the wheel colliding with constraints.

where ${\tilde{R}}_n$, ${\tilde{R}}_t$ —impulses described in the formulae (15); $\beta \in \left(0,1\right]$ —the coefficient of restitution, whose value for the elastic collision is 1.

4. Description of the Motion of the Wheel

The motion of the wheel will be described in the matrix form, using generalized coordinates. To this end, we introduce the following designations:

$X:=\left[\begin{array}{c}x\\ y\\ \phi \end{array}\right]$ —generalized coordinates vector,

$M:=\left[\begin{array}{ccc}m& 0& 0\\ 0& m& 0\\ 0& 0& J\end{array}\right]$ —wheel inertia matrix,

$Q:=\left[\begin{array}{c}0\\ -mg\\ 0\end{array}\right]$ —vector of the gravitational force acting on the wheel,

$F_n\in R^3$ —vector of the reaction force perpendicular to the boundary of the set $\Omega$,

$F_t\in R^3$ —vector of the reaction force tangent to the boundary of the set $\Omega$,

$X\in \Omega \subset R{}^3$ —relations describing constraining of the motion of the wheel.

If the position of the wheel is such that $X\in \text{Int}\Omega$, then the constraints do not act and the motion of the wheel is described as:

$M\text{}\ddot{X}=Q$, (17a)

and if the motion of the wheel is constrained, which may occur when $X\in \text{Fr}\Omega$, then the equation of motion should include a description of reaction forces

$M\text{}\ddot{X}=Q+F_n-F_t$. (17b)

In the next section of the paper we will define the method of determining the reaction forces.

4.1. Determination of Reaction Forces

In the descriptions of constraints and reaction forces presented further in the text we will use vector functions G and H, which to each vector X satisfying the relation $X\in \text{Fr}\Omega$ assign a pair of vectors, namely the perpendicular vector, $G\in R^3$ and the tangent vector, $H\in R^3$. With these vectors we can determine the velocities and acceleration of the point on the wheel in contact with the boundary of the set $\Omega$ (See Figure 3),

$v_n:=G^{\text{T}}\left(X\right)\stackrel{\dot{}}{X},v_t:=H^{\text{T}}\left(X\right)\stackrel{\dot{}}{X}$, (18a)

$a_n:=G^{\text{T}}\left(X\right)\ddot{X}+{\gamma }_n\left(X,\stackrel{\dot{}}{X}\right),a_t:=H^{\text{T}}\left(X\right)\ddot{X}+{\gamma }_t\left(X,\stackrel{\dot{}}{X}\right)$, (18b)

where ${\gamma }_n$ and ${\gamma }_t$ are scalar functions. Detailed forms of the vector functions, G and H, and scalar functions, ${\gamma }_n$ and ${\gamma }_t$, will be given further on in the paper.

Using the functions, G and ${\gamma }_n$ one can describe the consequence $D\Omega$, according to formula (6), determining the set of admissible velocities $\stackrel{\dot{}}{X}$,

$D\Omega \left(X\right):=\left\{\stackrel{\dot{}}{X}\in R^3:G^{\text{T}}\left(X\right)\stackrel{\dot{}}{X}\ge 0\right\}$, when $X\in \text{Fr}\Omega$ (19a)

and the set of admissible velocities according to the formula (8b),

$D^2\Omega \left(X,\stackrel{\dot{}}{X}\right):=\left\{\ddot{X}\in R^3:G^{\text{T}}\left(X\right)\ddot{X}+{\gamma }_n\left(X,\stackrel{\dot{}}{X}\right)\ge 0\right\}$,

when $X\in \text{Fr}\Omega$ and $G^{\text{T}}\left(X\right)\stackrel{\dot{}}{X}=0$ (19b)

If the position X and the velocity $\stackrel{\dot{}}{X}$ satisfy the relations specified in formula (19b), then the reaction force F occurs, whose two components $F_n$ and $F_t$ are present in the equation of motion (17b).

We will describe the generalized reaction forces $F_n,F_t\in R^3$ on the basis of the relations given in formulae (12). The normal reaction force is determined by the formula,

$F_n:=G\left(X\right)\lambda$, (20a)

if $\lambda \ge 0$ is the multiplier of the reaction satisfying relations resulting from the formula (12a), that is

$\lambda \ge 0,G^{\text{T}}\left(X\right)\ddot{X}+{\gamma }_n\left(X,\stackrel{\dot{}}{X}\right)\ge 0,\lambda \left(G^{\text{T}}\left(X\right)\ddot{X}+{\gamma }_n\left(X,\stackrel{\dot{}}{X}\right)\right)=0$ (20b)

If we take into account the description of relations from the formula (11), we will obtain an equivalent description,

$\lambda ={\left[\lambda -\rho \left(G^{\text{T}}\left(X\right)\ddot{X}+{\gamma }_n\left(X,\stackrel{\dot{}}{X}\right)\right)\right]}^{+}$, when $X\in \text{Fr}\Omega$ and $G^{\text{T}}\left(X\right)\stackrel{\dot{}}{X}=0$, (21)

if $\rho$ is any positive number.

Likewise, according to the formulae (12b) and (12c) we determine relations describing the tangential reaction $F_t\in R^3$,

$F_t=\mu H\left(X\right)\Vert G\Vert \lambda \tau$ (22a)

$\tau =\text{sign}{v}_t\equiv \text{sign}\left(H^{\text{T}}\left(X\right)\stackrel{\dot{}}{X}\right)$, when $v_t\ne 0$ (22b)

$\tau \in \{\begin{array}{ll}\left\{\text{sign}{a}_t\right\}\hfill & \text{for}{a}_t\ne 0\hfill \\ \left[-1,+1\right]\hfill & \text{for}{a}_t=0\hfill \end{array}$, when $v_t=0$, (22c)

if the vector $a_t$ has been described in the formulae (18b).

In the formula (20a), we have taken into consideration the fact that the load exerted by the wheel on the boundary of constraints, specified by the force from the formula (20a), is

$N:=\Vert F_n\Vert =\Vert G\Vert \lambda$ (23a)

and the resulting friction force limit (the Coulomb’s limit) is equal to

$T_{\text{o}}:=\mu N=\mu \Vert G\Vert \lambda$. (23b)

The relation between the multiplier $\tau$ and the acceleration $\ddot{X}$ given with the formula (22c) can be presented, analogically to the relation (21), in the form:

$\tau =\Pi \left(\tau +\rho \left(H^{\text{T}}\left(X\right)\ddot{X}+{\gamma }_t\left(X,\stackrel{\dot{}}{X}\right)\right)\right)$, when $\lambda >0$, (23c)

where $\rho$ is any positive number.

4.2. Wheel Acceleration

The equations of motion represent a set of relations based on which the wheel acceleration vector $\ddot{X}\left(t\right)\in R^3$ can be determined at any instant of time $t\in \left[0,t_{end}\right)$, provided that constraints of motion are satisfied, that is when $X\left(t\right)\in \Omega$ and $\stackrel{\dot{}}{X}\left(t\right)\in D\Omega \left(X\left(t\right)\right)$.

If the wheel is in contact with the boundary of constraint, that is when

$X\left(t\right)\in \text{Fr}\Omega$ and $\stackrel{\dot{}}{X}\left(t\right)\in \text{Fr}D\Omega \left(X\left(t\right)\right)$, (24)

then the constraint reaction forces considered in the previous section are acting on the wheel.

On this basis we formulate a set of relations from which the value of the acceleration vector can be unequivocally determined. The set includes:

- the equation of motion (17b), which includes the expressions (20a) and (22a)

$M\text{}\ddot{X}=Q+G\lambda -\mu H\Vert G\Vert \lambda \tau$ (25a)

- the relations specifying the multipliers according to the formulae (21) and (24)

$\lambda ={\left[\lambda -{\rho }_1\left(G^{\text{T}}\ddot{X}+{\gamma }_n\right)\right]}^{+},{\rho }_1>0$ (25b)

$\tau =\{\begin{array}{ll}\text{sign}\left(H^{\text{T}}\stackrel{\dot{}}{X}\right)\hfill & \text{for}{H}^{\text{T}}\stackrel{\dot{}}{X}\ne 0\hfill \\ \Pi \left(\tau +{\rho }_2\left(H^{\text{T}}\ddot{X}+{\gamma }_t\right)\right),{\rho }_2>0\hfill & \text{for}{H}^{\text{T}}\stackrel{\dot{}}{X}=0\hfill \end{array}$, (25c)

however, in the formulae (25) the arguments of the mappings G, H, ${\gamma }_n$ and ${\gamma }_t$ have not been specified.

Equation (25) determine at any moment in time $t\in \left[0,t_{end}\right)$, the value of the acceleration vector $\ddot{X}\left(t\right)$ and the values of multipliers $\lambda \left(t\right)$ and $\tau \left(t\right)$, which specify the values of the reaction forces.

The relations (25) comprise the mappings G, H, ${\gamma }_n$ and ${\gamma }_t$, the descriptions of which will be given in the next section.

4.3. Description of Additional Mappings

The mappings G, H, ${\gamma }_n$ and ${\gamma }_t$ given in the previous section of the paper serve as a description of the geometrical features of the boundary of constraints $\text{Fr}\Omega$, with the use of generalized coordinates. As mentioned before, the vector mapping G determines, in each point of the boundary $X\in \text{Fr}\Omega$, the vector $G\left(X\right)\in R^3$, which is perpendicular to to the boundary of constraints, and the mapping H determines the vector $H\left(X\right)$, which is tangent to it. The scalar mappings, ${\gamma }_n$ and ${\gamma }_t$ serve to determine the effect of the curvature of the boundary and of the velocity of the wheel on its acceleration.

The considered boundary of the set is presented in Figure 2. It is composed of four segments, including two singular points A and B (cf. Figure 7).

The description of the mappings under consideration depends on the curve segment and follows from the formulae (3) and (4). Hence, we obtain:

- Segment I, ${\text{Fr}}^{\text{I}}\Omega :=\left\{\left(x,y\right):x=r,y>r+h\right\}$

$G_{\text{I}}:=\left[\begin{array}{c}1\\ 0\\ 0\end{array}\right],H_{\text{I}}:=\left[\begin{array}{c}0\\ 1\\ -r\end{array}\right],{\gamma }_n^{\text{I}}={\gamma }_t^{\text{I}}=0$ (26a)

- Segment II, ${\text{Fr}}^{\text{II}}\Omega :=\left\{\left(x,y\right):r<x\le l_1,y=r+h\right\}$

$G_{\text{II}}:=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right],H_{\text{II}}:=\left[\begin{array}{c}1\\ 0\\ r\end{array}\right],{\gamma }_n^{\text{II}}={\gamma }_t^{\text{II}}=0$ (26b)

- Segment III, ${\text{Fr}}^{\text{III}}\Omega :=\left\{\left(x,y\right):l_1<x<l_1+a,y=f\left(x\right)\right\}$

$G_{\text{III}}:=\left[\begin{array}{c}-f^{\prime }\left(x\right)\\ 1\\ 0\end{array}\right],H_{\text{III}}:=\left[\begin{array}{c}1\\ f^{\prime }\left(x\right)\\ r\end{array}\right],{\gamma }_n^{\text{III}}:=f^{″}\left(x\right){\stackrel{\dot{}}{x}}^2,{\gamma }_t^{\text{III}}:=f^{″}\left(x\right)\stackrel{\dot{}}{x}\stackrel{\dot{}}{y}$ (26c)

- Segment IV, ${\text{Fr}}^{\text{IV}}\Omega :=\left\{\left(x,y\right):x>l_1+a,y=r\right\}$

$G_{\text{IV}}:=\left[\begin{array}{c}0\\ 1\\ 0\end{array}\right],H_{\text{IV}}:=\left[\begin{array}{c}1\\ 0\\ r\end{array}\right],{\gamma }_n^{\text{IV}}={\gamma }_t^{\text{IV}}=0$ (26d)

In the formulae (26) we have used designations marked in Figure 2 and associated with the function f, described in the formula (4c). Thus, we obtain:

$\begin{array}{l}{f}^{\prime }\left(x\right):=-\frac{x-l_1}{\sqrt{r^2-{\left(x-l_1\right)}^2}}\hfill \\ f^{″}\left(x\right):=-\frac{r^2}{{\left[r^2-{\left(x-l_1\right)}^2\right]}^{3/2}}\hfill \end{array}\}\text{for}{l}_1<x\le l_1+a$ (26e)

The singular positions of the wheel are specified by the wheel centre point coordinates, i.e.

Figure 7. The boundary of the set $\text{Fr}\Omega$ with marked segments and singular points.

$\text{point}A:x=r,y=r+h$ (27a)

$\text{point}B:x=l_1+a,y=r$. (27b)

In these positions the wheel comes into contact with the constraints at two points (Figure 8).Thus, the resulting descriptions of the sets of admissible velocities are as follows:

$D\Omega \left(X_A\right):=\left\{\stackrel{\dot{}}{X}\in R^3:G_{\text{I}}^{\text{T}}\stackrel{\dot{}}{X}\ge 0,G_{\text{II}}^{\text{T}}\stackrel{\dot{}}{X}\ge 0\right\}$ (28a)

$D\Omega \left(X_B\right):=\left\{\stackrel{\dot{}}{X}\in R^3:G_{\text{II}}^{\text{T}}\stackrel{\dot{}}{X}\ge 0,G_{\text{III}}^{\text{T}}\stackrel{\dot{}}{X}\ge 0\right\}$, (28b)

where designations from formulae (26) have been used.

As mentioned before, a constraint to the wheel’s acceleration and a reaction force arise when the velocity vector is situated at the boundary of the set $D\Omega$, which has been described in the formulae (19) and (20).

On this basis, we specify the relations for the normal component of the reaction force in the position A (27a)

$F_n^A:=\{\begin{array}{lll}{G}_{\text{I}}{\lambda }_{\text{I}}\hfill & \text{if}{G}_{\text{I}}^{\text{T}}\stackrel{\dot{}}{X}=0,G_{\text{II}}^{\text{T}}\stackrel{\dot{}}{X}>0\hfill & \text{forpoint}{A}_1\hfill \\ G_{\text{II}}{\lambda }_{\text{II}}\hfill & \text{if}{G}_{\text{I}}^{\text{T}}\stackrel{\dot{}}{X}>0,G_{\text{II}}^{\text{T}}\stackrel{\dot{}}{X}=0\hfill & \text{forpoint}{A}_2\hfill \\ G_{\text{I}}{\lambda }_{\text{I}}+G_{\text{II}}{\lambda }_{\text{II}}\hfill & \text{if}{G}_{\text{I}}^{\text{T}}\stackrel{\dot{}}{X}=0,G_{\text{II}}^{\text{T}}\stackrel{\dot{}}{X}=0\hfill & \text{forbothpoints}\hfill \end{array}$. (29)

Values of the multipliers ${\lambda }_{\text{I}}$ and ${\lambda }_{\text{II}}$ are connected with the vector of acceleration $\ddot{X}$ by the relation described in the formula (25b).

The force of reaction $F_n^B$, acting on the wheel in position B (27b) is calculated analogically to the formula (29), i.e.

$F_n^B:=\{\begin{array}{lll}{G}_{\text{III}}{\lambda }_{\text{III}}\hfill & \text{if}{G}_{\text{III}}^{\text{T}}\stackrel{\dot{}}{X}=0,G_{\text{IV}}^{\text{T}}\stackrel{\dot{}}{X}>0\hfill & \text{forpoint}{B}_1\hfill \\ G_{\text{IV}}{\lambda }_{\text{IV}}\hfill & \text{if}{G}_{\text{III}}^{\text{T}}\stackrel{\dot{}}{X}>0,G_{\text{IV}}^{\text{T}}\stackrel{\dot{}}{X}=0\hfill & \text{forpoint}{B}_2\hfill \\ G_{\text{III}}{\lambda }_{\text{III}}+G_{\text{IV}}{\lambda }_{\text{IV}}\hfill & \text{if}{G}_{\text{III}}^{\text{T}}\stackrel{\dot{}}{X}=0,G_{\text{IV}}^{\text{T}}\stackrel{\dot{}}{X}=0\hfill & \text{forbothpoints}\hfill \end{array}$ (30)

The form of the description of the relations between the multipliers $\lambda$ and the vector of acceleration $\ddot{X}$ has been specified by the formula (25b).

The tangential component of the reaction force $F_t$ will be specified based on

Figure 8. The singular positions in which the wheel makes contact with the constraints at two points.

the relation (22), but after taking into account the proper indices of the vectors G, H and of the multipliers $\lambda$, $\tau$.

5. Description of the Collision

In the situation specified by the relation (9), the velocity of the wheel is inadmissible due to constraints, which results in a collision. As a result, there occurs an impulse reaction of the constraints, adjusting the wheel’s velocity to the limitations. In section 3 we have presented a local description of the collision model, in the form of a relation specifying the association between the impulse reaction and the abrupt change in the velocity of the wheel. We will now present the relations described in that section using the generalized coordinates.

Description of a plastic collision includes also the momentum balance equation and the relations specifying the impulse reaction forces. The general form of this description is analogical to the equations of motion (17) and (25), i.e.

$M\text{}{\stackrel{\dot{}}{X}}^{+}-M\text{}{\stackrel{\dot{}}{X}}^{-}={\tilde{F}}_n-{\tilde{F}}_t$ (31)

${\tilde{F}}_n=G\tilde{\lambda },{\tilde{F}}_t=\mu \Vert G\Vert H\tilde{\lambda }\tau$ (32)

and the description of the relations from the formulae (15) now has a form analogical to (25b) and (25c), i.e.

$\tilde{\lambda }={\left[\tilde{\lambda }-{\rho }_1{G}^{\text{T}}{\stackrel{\dot{}}{X}}^{+}\right]}^{+}$ (33a)

$\tau =\Pi \left(\tau +{\rho }_2{H}^{\text{T}}{\stackrel{\dot{}}{X}}^{+}\right)$, (33b)

where the force impulses acting on the wheel have been marked with the index “~”.

When the collision of the wheel with the constraints occurs in the singular position of the wheel specified by point A or B (according to the formulae (27)), then two reaction forces are acting on the wheel, i.e.

${\tilde{F}}_n=G_{\text{I}}{\tilde{\lambda }}_{\text{I}}+G_{\text{II}}{\tilde{\lambda }}_{\text{II}}$ for point A, (34a)

${\tilde{F}}_n=G_{\text{III}}{\tilde{\lambda }}_{\text{III}}+G_{\text{IV}}{\tilde{\lambda }}_{\text{IV}}$ for point B. (34b)

The relations specifying the multipliers, $\tilde{\lambda }$ and $\tau$, from the formulae (34), have the same form as in the formulae (33) but after taking into account the suitable indices.

As mentioned before, the relations from the formulae (32) and (33) specify an abrupt change in the velocity of the wheel $\Delta V:=V^{+}-V^{-}$, and the reaction force impulses, ${\tilde{F}}_n$ and ${\tilde{F}}_t$, arising during the plastic collision.

If an elastic or elasto-plastic collision is considered, then—according to the Newton’s and Poisson’s hypothesis specified by the formulae (16)—we need to take into account an additional collision phase (restitution) described by the equations

$M\text{}{\stackrel{\dot{}}{X}}^{\ast }-M\text{}{\stackrel{\dot{}}{X}}^{+}={\tilde{F}}_n^{\ast }-{\tilde{F}}_t^{\ast }$, (35a)

where ${\stackrel{\dot{}}{X}}^{\ast }\in R^3$ —the velocity vector after the elastic or elasto-plastic collision; $F_n^{\ast }$, $F_t^{\ast }$ —the reaction force impulses generated in this phase, specified by the relations:

$F_n^{\ast }=G{\lambda }^{\ast }$, ${\lambda }^{\ast }=\beta \tilde{\lambda }$, (35b)

$F_t^{\ast }=\mu \Vert G\Vert H{\lambda }^{\ast }{\tau }^{\ast }$, ${\tau }^{\ast }=\Pi \left({\tau }^{\ast }+{\rho }_2{H}^{\text{T}}{\stackrel{\dot{}}{X}}^{\ast }\right)$. (35c)

In the above relations, the value of the multiplier ${\lambda }^{\ast }$ is known, which follows from the formula (35c), where $\beta \in \left(0,1\right]$ is the coefficient of restitution, the value of which for the elastic collision is 1.

Relations (35) determine the velocity vector after the collision ${\stackrel{\dot{}}{X}}^{\ast }$ and the multiplier ${\tau }^{\ast }$, specifying the effect of friction on this collision phase.

6. Wheel Motion Simulation

6.1. Description of the Method of Solving the Initial Problem

The non-smooth problem of a mechanical system motion was closely analyzed in the monograph [1], where perfect unilateral constraints were considered. It follows from these considerations that the motion of such system is described by an absolutely continuous function. A staged method of determining such a function was also specified in this work. We will use this method to determine the function describing the motion of the wheel.

The descriptions of the relations considered in sections 4 and 5, specifying the motion of the wheel will be now presented in a synthetic form by means of two equations,

$\ddot{X}=D\left(X,\stackrel{\dot{}}{X}\right),\text{if}X\in \Omega ,\stackrel{\dot{}}{X}\in D\Omega \left(X\right)$ (36)

${\stackrel{\dot{}}{X}}^{+}=D\left({\stackrel{\dot{}}{X}}^{-},X\right),\text{if}X\in \text{Fr}\Omega ,{\stackrel{\dot{}}{X}}^{-}\notin D\Omega \left(X\right)$, (37)

where $X\in R^3$ —the vector of the generalized coordinates, $D:R^6\to R^3$ —the mapping determining the vector of acceleration, $D:R^3\to R^3$ —the mapping determining an abrupt change in velocity in the moments of collisions.

The mapping $A$ is specified in the implicit form by means of the relation (25). In this case, determining the value of the mapping $A$ comes down to determining the solution of the relation (25), i.e. to calculating the values of multipliers $\lambda$ and $\tau$, and of the vector $\ddot{X}$.

The mapping $A$ is also implicit and its values are determined by the solution to relations (32)-(35), specified by the multipliers $\tilde{\lambda }$ and $\tau$ and the vector ${\stackrel{\dot{}}{X}}^{+}$.

At every stage of the above-mentioned method the solution of the initial problem specified by Equation (36) is determined. The continuous and differentiable function X so established, describes the motion of the wheel up to the moment when the collision specified by the conditions from the formula (37) occurs. Then an abrupt change in the wheel’s velocity described by the mapping $A$ takes place. In this situation the considered stage of the solution comes to an end.

This staged method of finding the initial problem solution has served as a basis for developing calculation software for wheel motion simulation. The program includes algorithms for determining the values of implicit mappings $A$ and $A$, in all segments of the set $\text{Fr}\Omega$ described by the formulae (26) and at all singular points (27) designated in Figure 7. It needs to be added, that in order to determine the values of the implicit mapping $A$, it is necessary to solve a system of algebraic equations containing five unknowns: $\ddot{X}\in R^3$, $\lambda$ and $\tau$. The values of the implicit mapping $A$ are determined in the same way.

6.2. Detailed form of the Equations of Motion and Reaction

To start with, we will consider the detailed form of equations of motion and reaction, when the wheel makes contact with the boundary of constraints in sector IV (See Figure 7). In Figure 3 the wheel with index 6 illustrates such a position. In the considered situation, the conditions specified in the formula (24), that is $x>l_1+a$, $y=r$ and $\stackrel{\dot{}}{y}=0$ are satisfied, where a, $l_1$ —dimensions described in the formula (4b).

The equations of motion (25) in expanded form, and after taking into account the formulae (26d), for sector IV, have the form:

$m\ddot{x}=-\mu \lambda \tau$ (38)

$m\ddot{y}=-mg+\lambda$ (39)

$J\ddot{\phi }=-\mu \lambda r\tau$ (40)

$\lambda ={\left[\lambda -\rho \ddot{y}\right]}^{+},\rho >0$ (41)

$\tau =\{\begin{array}{ll}\text{sign}\left(\stackrel{\dot{}}{x}+r\stackrel{\dot{}}{\phi }\right)\hfill & \text{if}\stackrel{\dot{}}{x}+r\stackrel{\dot{}}{\phi }\ne 0\hfill \\ \Pi \left(\tau +\rho \left(\ddot{x}+r\ddot{\phi }\right)\right),\rho >0\hfill & \text{if}\stackrel{\dot{}}{x}+r\stackrel{\dot{}}{\phi }=0\hfill \end{array}$. (42)

The above system of equations describes the accelerations of the wheel ( $\ddot{x}$, $\ddot{y}$, $\ddot{\phi }$ ), and the multipliers of reactions ( $\lambda$, $\tau$ ). This is an implicit form of the mapping $A$, when the wheel is in contact with the boundary of sector IV. Determining the solution to the above system of equations is relatively simple:

- if the wheel is slipping, i.e. $\stackrel{\dot{}}{x}+r\stackrel{\dot{}}{\phi }\ne 0$, then the value $\tau$ is determined directly from Equation (42₁);

- if the wheel is rolling without slip, then from Equations (42₂) and (38) and (39) we obtain that $\tau =0$ ;

- from the relatively simple system of Equations (39) and (41) we obtain $\lambda =mg$ ;

- having determined the values, $\lambda$ and $\tau$, from Equations (38), (39) and (40), we can determine acceleration values, and in particular we have $\ddot{y}=0$.

The above-described Equations (38)-(42) apply to motion of the wheel moving (rolling) on the boundary. Now, we will deal with describing a collision of the wheel with the boundary in sector IV.A general description of the collision has been given in the formulae (32), (33) and (35). In the considered situation the expanded description of the collision, assuming the plastic collision hypothesis, has the form:

$m{\stackrel{\dot{}}{x}}^{+}-m{\stackrel{\dot{}}{x}}^{-}=-\mu \tilde{\lambda }\tau$ (43)

$m{\stackrel{\dot{}}{y}}^{+}-m{\stackrel{\dot{}}{y}}^{-}=\tilde{\lambda }$ (44)

$J{\stackrel{\dot{}}{\phi }}^{+}-J{\stackrel{\dot{}}{\phi }}^{-}=-\mu \tilde{\lambda }r\tau$ (45)

$\tilde{\lambda }={\left[\tilde{\lambda }-\rho {\stackrel{\dot{}}{x}}^{+}\right]}^{+},\rho >0$ (46)

$\tau =\Pi \left(\tau +\rho v_t^{+}\right),\rho >0$ (47)

$v_t^{+}:={\stackrel{\dot{}}{x}}^{+}+r{\stackrel{\dot{}}{\phi }}^{+}$. (48)

This is a system of equations in regard to ${\stackrel{\dot{}}{x}}^{+}$, ${\stackrel{\dot{}}{y}}^{+}$, $\tilde{\lambda }$, $\tau$. The solution is determined in the following steps:

- from Equations (44) and (46) we can directly determine $\tilde{\lambda }=-m{\stackrel{\dot{}}{y}}^{-}$ ;

- according to the formula (48) and Equations (43) and (44) we obtain the relationship,

$v_t^{+}:=v_t^{-}-\left(\frac{1}{m}+\frac{r^2}{J}\right)\mu \tilde{\lambda }\tau ;v_t^{-}:={\stackrel{\dot{}}{x}}^{-}+r{\stackrel{\dot{}}{\phi }}^{-}$, (49)

and next, from Equation (47) we determine the values, $\tau$ and $v_t^{+}$ ;

- having determined the values, $\tilde{\lambda }$ and $\tau$ we calculate the post-collision velocities, ${\stackrel{\dot{}}{x}}^{+}$ and ${\stackrel{\dot{}}{y}}^{+}$.

The solution determined specifies the effect of the plastic collision. When considering an elastic or elasto-plastic collision it is necessary to take into account the change in velocity in the second phase of the collision, described in the formulae (35). In the example under consideration, this description, with respect to the coefficient of restitution $\beta \in \left(0,1\right]$, has the form:

${\tilde{\lambda }}^{\ast }=\beta \tilde{\lambda }$, (50)

$m{\stackrel{\dot{}}{x}}^{\ast }=m{\stackrel{\dot{}}{x}}^{+}-\mu {\tilde{\lambda }}^{\ast }{\tau }^{\ast }$, (51)

$m{\stackrel{\dot{}}{y}}^{\ast }=m{\stackrel{\dot{}}{y}}^{+}+{\tilde{\lambda }}^{\ast }$, (52)

$J{\stackrel{\dot{}}{\phi }}^{\ast }=J{\stackrel{\dot{}}{\phi }}^{+}-\mu {\tilde{\lambda }}^{\ast }r{\tau }^{\ast }$ (53)

${\tau }^{\ast }=\Pi \left({\tau }^{\ast }+\rho v_t^{\ast }\right),\rho >0$, (54)

$v_t^{\ast }:={\stackrel{\dot{}}{x}}^{\ast }+r{\stackrel{\dot{}}{\phi }}^{\ast }=v_t^{+}-\left(\frac{1}{m}+\frac{r^2}{J}\right)\mu {\tilde{\lambda }}^{\ast }{\tau }^{\ast }$. (55)

The effect of the collision is described by the velocities ${\stackrel{\dot{}}{x}}^{\ast }$, ${\stackrel{\dot{}}{y}}^{\ast }$, ${\stackrel{\dot{}}{\phi }}^{\ast }$ and multipliers ${\tilde{\lambda }}^{\ast }$ and ${\tau }^{\ast }$. The value ${\tilde{\lambda }}^{\ast }$ is specified by the formula (50), while the values ${\tau }^{\ast }$ and $v_t^{+}$ are obtained directly from the formulae (54) and (55).

During the motion of the wheel in sector III (cf. Figure 6), the wheel makes contact only with the corner of the curb, which is illustrated in Figure 3 (the wheel with index 4). In this case the centre point of the wheel takes the position specified by the conditions (see formula (26c)),

$l_1\le x<l_1+a$, $y=f\left(x\right)$ (56a)

and the velocity of the wheel is such, that

$v_n:=\stackrel{\dot{}}{y}-f^{\prime }\left(x\right)\stackrel{\dot{}}{x}=0$. (56b)

The condition (56b) indicates, that the normal component of the velocity of the point on the wheel (Figure 3 and formula (10)) making contact with the constraints equals zero, i.e. the wheel does not get separated from the curb.

In the formulae (26c) we included definitions of the vectors G and H and expressions ${\gamma }_n$ and ${\gamma }_t$. To simplify further descriptions we will omit indices “III” present next to the enumerated designations.

For the wheel status so specified (formulae (56a) and (56b)), we formulate an expanded description of the wheel’s acceleration and constraint reactions according to the formulae (25):

$m\ddot{x}=-f^{\prime }\left(x\right)\lambda -\mu \Vert G\Vert \lambda \tau$, (57)

$m\ddot{y}=-mg+\lambda -\mu \Vert G\Vert f^{\prime }\left(x\right)\lambda \tau$, (58)

$J\ddot{\phi }=-\mu \Vert G\Vert \lambda r\tau$, (59)

$\lambda ={\left[\lambda -\rho a_n\right]}^{+},\rho >0$, (60)

$\tau =\{\begin{array}{ll}\text{sign}{v}_t\hfill & \text{if}{v}_t\ne 0\hfill \\ \Pi \left(\tau +\rho a_t\right),\rho >0\hfill & \text{if}{v}_t=0\hfill \end{array}$, (61)

where the following designations have been introduced:

$G\left(x\right):={\left[\begin{array}{ccc}-f^{\prime }\left(x\right)& 1& 0\end{array}\right]}^{\text{T}}$, (according to the formula (26c)) (62)

$\Vert G\Vert \lambda$ —the contact force exerted by the wheel on the curb,

$v_t$ —the tangential component of the velocity of the point on the wheel in contact with the curb,

$v_t:=\stackrel{\dot{}}{x}+f^{\prime }\left(x\right)\stackrel{\dot{}}{y}+r\stackrel{\dot{}}{\phi }$ (63)

$a_n$ —the normal component of the acceleration of the point on the wheel in contact with the curb

$a_n:=\ddot{y}-f^{\prime }\left(x\right)\ddot{x}+f^{″}\left(x\right){\stackrel{\dot{}}{x}}^2$, (64a)

from the condition of constraints (11) it follows that the value of this acceleration may be greater than or equal to zero. If in the formula (64a) we include Equations (57) and (58), we will obtain:

$a_n:=\frac{1}{m}\left({\Vert G\Vert }^2\lambda +\left(\gamma -g\right)\right),\gamma :=f^{″}\left(x\right){\stackrel{\dot{}}{x}}^2$ (64b)

$a_t$ —the tangential component of the acceleration of the point on the wheel n contact with the curb,

$a_t:=\ddot{x}+f^{\prime }\left(x\right)\ddot{y}+r\ddot{\phi }+f^{″}\left(x\right)\stackrel{\dot{}}{x}\stackrel{\dot{}}{y}$, (65a)

if in the above equation we include Equations (57), (58) and (59), and make suitable transformations, we will obtain:

$a_t:=-\left(\frac{1}{m}{\Vert G\Vert }^2+\frac{r^2}{J}\right)\mu \Vert G\Vert \lambda \tau +f^{\prime }\left(x\right)g+f^{″}\left(x\right)\stackrel{\dot{}}{x}\stackrel{\dot{}}{y}$. (65b)

To determine the value of multiplier $\lambda$ we need to substitute the expression $a_n$ into Equation (60) according to the formula (64b); from this we can directly calculate the value $\lambda$.

When the value of the velocity of the wheel’s slip on the curb $v_t$ (according to the formula (63)) is different from zero, then from the formula (61₁) we obtain the value of the multiplier $\tau =\text{sign}{v}_t$. If $v_t=0$, then we should substitute the expression $a_t$ into Equation (61₂) according to the formula (65b); in this way we obtain the value of the multiplier $\tau$.

Having determined the values of the multipliers $\lambda$ and $\tau$ we calculate the values of the accelerations $\ddot{x}$, $\ddot{y}$ and $\ddot{\phi }$ according to the formulae (57), (58) and (59).

The collision of the wheel with the curb occurs when the velocity $v_n$ according to the formula (56b) is negative, i.e. $v_n<0$ (See Figure 3(a)). The effect of this phenomenon, on the assumption of the plastic collision hypothesis, is described by the equations:

$m{\stackrel{\dot{}}{x}}^{+}-m{\stackrel{\dot{}}{x}}^{-}=-f^{\prime }\left(x\right)\tilde{\lambda }-\mu \Vert G\Vert \tilde{\lambda }\tau$ (66)

$m{\stackrel{\dot{}}{y}}^{+}-m{\stackrel{\dot{}}{y}}^{-}=\tilde{\lambda }-\mu f^{\prime }\left(x\right)\Vert G\Vert \tilde{\lambda }\tau$, (67)

$J{\stackrel{\dot{}}{\phi }}^{+}-J{\stackrel{\dot{}}{\phi }}^{-}=-\mu \Vert G\Vert \tilde{\lambda }r\tau$, (68)

$\tilde{\lambda }={\left[\tilde{\lambda }-\rho v_n^{+}\right]}^{+},\rho >0$, (69)

$\tau =\Pi \left(\tau +\rho v_t^{+}\right),\rho >0$. (70)

Solving these equations in regard to ${\stackrel{\dot{}}{x}}^{+}$, ${\stackrel{\dot{}}{y}}^{+}$, ${\stackrel{\dot{}}{\phi }}^{+}$, $\tilde{\lambda }$, $\tau$ determines the velocity of the wheel after the plastic collision and the reaction force impulses. To determine the solution, we establish the expressions specifying $v_n^{+}$ and $v_t^{+}$, depending on the multipliers, $\tilde{\lambda }$ and $\tau$.

$v_n^{+}:={\stackrel{\dot{}}{y}}^{+}-f^{\prime }\left(x\right){\stackrel{\dot{}}{x}}^{+}=v_n^{-}+\frac{1}{m}{\Vert G\Vert }^2\tilde{\lambda }$ (71)

$v_t^{+}:={\stackrel{\dot{}}{x}}^{+}-f^{\prime }\left(x\right){\stackrel{\dot{}}{y}}^{+}+r{\stackrel{\dot{}}{\phi }}^{+}=v_t^{-}+\left(\frac{1}{m}{\Vert G\Vert }^2+\frac{r^2}{J}\right)\mu \Vert G\Vert \tilde{\lambda }\tau$ (72)

Having included the formula (71) in Equation (69), we obtain the value of multiplier $\tilde{\lambda }$. Next, from Equation (70) and from the relationship (72) we determine the values, $\tau$ and $v_t^{+}$. Having determined the multipliers we calculate the values of velocities looked for after the plastic collision, ${\stackrel{\dot{}}{x}}^{+}$, ${\stackrel{\dot{}}{y}}^{+}$ and ${\stackrel{\dot{}}{\phi }}^{+}$.

If we consider an elastic or elasto-plastic collision, we should consider a description of the second phase of collision (35), analogical to the formulae (50)-(55),

${\tilde{\lambda }}^{\ast }=\beta \tilde{\lambda }$, $\beta \in \left(0,1\right]$ (73)

$m{\stackrel{\dot{}}{x}}^{\ast }=m{\stackrel{\dot{}}{x}}^{+}-f^{\prime }\left(x\right){\tilde{\lambda }}^{\ast }-\mu \Vert G\Vert {\tilde{\lambda }}^{\ast }{\tau }^{\ast }$ (74)

$m{\stackrel{\dot{}}{y}}^{\ast }=m{\stackrel{\dot{}}{y}}^{+}-\mu f^{\prime }\left(x\right)\Vert G\Vert {\tilde{\lambda }}^{\ast }{\tau }^{\ast }$, (75)

$J{\stackrel{\dot{}}{\phi }}^{\ast }=J{\stackrel{\dot{}}{\phi }}^{+}-\mu \Vert G\Vert {\tilde{\lambda }}^{\ast }r{\tau }^{\ast }$, (76)

${\tau }^{\ast }=\Pi \left({\tau }^{\ast }+\rho v_t^{\ast }\right),\rho >0$, (77)

$v_t^{\ast }:=v_t^{+}+\left(\frac{1}{m}{\Vert G\Vert }^2+\frac{r^2}{J}\right)\mu \Vert G\Vert {\tilde{\lambda }}^{\ast }{\tau }^{\ast }$. (78)

Based on Equations (77) and (78) we can directly determine the value of the multiplier ${\tau }^{\ast }$, and subsequently the values of the wheel’s velocity.

If the wheel assumes the singular position specified by point B (Figure 8), then the description of the action on the wheel becomes more complex, as the wheel makes contact with constraints at two points B₁ and B₂ marked in Figure 8. The position of the wheel is determined by the coordinates of the wheel’s centre point, $x=l_1+a$, $y=r$, while the description of the set of the admissible velocities of the wheel at this position is as follows:

$D\Omega \left(X_B\right):=\left\{\stackrel{\dot{}}{x},\stackrel{\dot{}}{y}\in R^1:-{f^{\prime }}_B\stackrel{\dot{}}{x}+\stackrel{\dot{}}{y}\ge 0,\stackrel{\dot{}}{y}\ge 0\right\}$, (79)

where ${f^{\prime }}_B$ —the value of the derivative of the function f for $x=l_1+a$

${f^{\prime }}_B:=-\frac{a}{\sqrt{r{}^2-a^2}}$, $a<r$ (80)

Figure 9 shows a schematic graph of the set $D\Omega$.

It is clear from an analysis of Figure 9, that in the considered position the constraints are acting, that is they are limiting accelerations, only when the velocity vector belongs to the boundary of the set $D\Omega$, which is written in the formula (8a). Figure 9 shows three examples of velocity vectors $V:={\left[\stackrel{\dot{}}{x},\stackrel{\dot{}}{y}\right]}^{\text{T}}$ situated on the boundary. For the vectors $V_1$ and $V_2$, only one constraint is acting, while for the vector $V_0$ both constraints are acting.

The action of individual constraints has been described before, and now we will deal with the situation specified by velocity $V_0$ that is, situation when $\stackrel{\dot{}}{y}-{f^{\prime }}_B\stackrel{\dot{}}{x}=0$ and $\stackrel{\dot{}}{y}=0$. It follows from these conditions, that the translational

Figure 9. Graph of the set $D\Omega$ in the singular point B (cf. Figure 7).

velocity equals zero, but we assume that the wheel is rotating, i.e. $\stackrel{\dot{}}{\phi }\ne 0$. In this situation the wheel is slipping on the constraints.

In the state described above the acceleration of the wheel is constrained, and the set of admissible accelerations has the form

$D^2\Omega \left(X_B,V_B\right):=\left\{\ddot{x},\ddot{y}\in R^1:-{f^{\prime }}_B\ddot{x}+\ddot{y}\ge 0,\ddot{y}\ge 0\right\}$. (81)

The equations specifying the acceleration and reaction forces now have the form:

$m\ddot{x}=-{f^{\prime }}_B{\lambda }_1-\mu \Vert G\Vert {\lambda }_1{\tau }_1-\mu {\lambda }_2{\tau }_2$ (82)

$m\ddot{y}=-mg+{\lambda }_1-\mu \Vert G\Vert {f^{\prime }}_B{\lambda }_1{\tau }_1+{\lambda }_2$, (83)

, (84)

${\lambda }_1={\left[{\lambda }_1-\rho a_1^n\right]}^{+},\rho >0$, (85)

${\lambda }_2={\left[{\lambda }_2-\rho a_2^n\right]}^{+},\rho >0$, (86)

where the reaction multipliers acting on the wheel at points $B_i$ showed in Figure 7 are marked as ${\lambda }_i$, ${\tau }_i$, for $i\in \left\{1,2\right\}$ ; and $a_1^n$, $a_2^n$ —normal components of the acceleration of the wheel at points $B_1$ and $B_2$.

$a_1^n:=\ddot{y}-{f^{\prime }}_B\ddot{x}$, $a_2^n:=\ddot{y}$ (87)

Since we have assumed that $\stackrel{\dot{}}{\phi }\ne 0$, then

${\tau }_1={\tau }_2=\text{sign}\stackrel{\dot{}}{\phi }$ (88)

On the basis of the above equations we should determine the values $\ddot{x}$, $\ddot{y}$, $\ddot{\phi }$, and ${\lambda }_1$ and ${\lambda }_2$. To this end, from the formulae (87) and Equations (82), (83) and (84) we establish the expressions:

$a_1^n=-g+\frac{1}{m}\left[{\Vert G\Vert }^2{\lambda }_1+\left(\mu {f^{\prime }}_B{\tau }_2+1\right){\lambda }_2\right]$. (89)

$a_2^n=-g+\frac{1}{m}\left[\left(1-\mu {f^{\prime }}_B\Vert G\Vert {\tau }_1\right){\lambda }_1+{\lambda }_2\right]$. (90)

Upon substituting the expression (89) into the formula (85) and assuming $\rho :=m/{\Vert G\Vert }^2$, we obtain

${\lambda }_1={\Vert G\Vert }^{-2}{\left[mg-{\gamma }_2{\lambda }_2\right]}^{+}$, if ${\gamma }_2:=\mu {f^{\prime }}_B{\tau }_2+1$. (91)

Likewise, from Equations (90) and (86), when $\rho :=m$, we obtain

${\lambda }_2={\left[mg-{\gamma }_1{\lambda }_1\right]}^{+}$, if ${\gamma }_1:=1-\mu {f^{\prime }}_B\Vert G\Vert {\tau }_1$. (92)

The solution of Equations (91) and (92) determines the values of multipliers ${\lambda }_1$ and ${\lambda }_2$, while the values of accelerations $\ddot{x}$, $\ddot{y}$, $\ddot{\phi }$ are calculated from the formulae (82), (83) and (84) after taking consideration (88).

It follows from an analysis of Equations (91) and (92) that, if

$\stackrel{\dot{}}{\phi }<0$, i.e. ${\tau }_1={\tau }_2=-1$, then (93a)

${\lambda }_1=0$, ${\lambda }_2=mg$, (93b)

that is an obvious result is obtained—the wheel moves away from the curb, that is $\ddot{x}=-\mu g$, $\ddot{y}=0$, $\ddot{\phi }=-\frac{1}{J}\mu mgr$.

If $\stackrel{\dot{}}{\phi }>0$, i.e. ${\tau }_1={\tau }_2=1$, (94a)

then two solutions are possible:

${\lambda }_1>0$, ${\lambda }_2>0$, when ${\Vert G\Vert }^2>{\gamma }_1$, (94b)

${\lambda }_1=\frac{mg}{\Vert G\Vert {}^2}$, ${\lambda }_2=0$, when ${\Vert G\Vert }^2\le {\gamma }_1$. (94c)

In the case specified by the formula (94b), the centre point of the wheel is not changing position $\ddot{x}=\ddot{y}=0$, but the wheel is rotating $\stackrel{\dot{}}{\phi }>0$, $\ddot{\phi }<0$ (spinning); while in the second case (94c) the wheel is climbing the curb, because $\ddot{x}<0$, $\ddot{y}>0$.

In the considered position of the wheel, a collision with the curb occurs, when the velocities ${\stackrel{\dot{}}{x}}^{-}$ and ${\stackrel{\dot{}}{y}}^{-}$ fail to satisfy the conditions given in the formula (79), that is they do not belong to the set $D\Omega \left(X_B\right)$ whose graph is shown in Figure 8. The resulting plastic collision is described by the following relations, given in Section 5, in the formulae (32)-(34):

$m{\stackrel{\dot{}}{x}}^{+}=m{\stackrel{\dot{}}{x}}^{-}-{f^{\prime }}_B{\tilde{\lambda }}_1-\mu \Vert G\Vert {\tilde{\lambda }}_1{\tau }_1-\mu {\tilde{\lambda }}_2{\tau }_2$, (95)

$m{\stackrel{\dot{}}{y}}^{+}=m{\stackrel{\dot{}}{y}}^{-}+{\tilde{\lambda }}_1-\mu {f^{\prime }}_B\Vert G\Vert {\tilde{\lambda }}_1{\tau }_1+{\tilde{\lambda }}_2$, (96)

$J{\stackrel{\dot{}}{\phi }}^{+}=J{\stackrel{\dot{}}{\phi }}^{-}-\mu \Vert G\Vert r{\tilde{\lambda }}_1{\tau }_1-\mu r{\tilde{\lambda }}_2{\tau }_2$, (97)

${\tilde{\lambda }}_1={\left[{\tilde{\lambda }}_1-\rho v_{n1}^{+}\right]}^{+},\rho >0$, (98)

${\tilde{\lambda }}_2={\left[{\tilde{\lambda }}_2-\rho v_{n2}^{+}\right]}^{+},\rho >0$, (99)

${\tau }_1=\Pi \left({\tau }_1+\rho v_{t1}^{+}\right),\rho >0$, (100)

${\tau }_2=\Pi \left({\tau }_2+\rho v_{t2}^{+}\right),\rho >0$, (101)

where the following designations have been introduced:

$v_{n1}^{-}:={\stackrel{\dot{}}{x}}^{-}-{f^{\prime }}_B{\stackrel{\dot{}}{y}}^{-}$, $v_{n1}^{+}:={\stackrel{\dot{}}{x}}^{+}-{f^{\prime }}_B{\stackrel{\dot{}}{y}}^{+}$, (102)

$v_{n2}^{-}:={\stackrel{\dot{}}{y}}^{-}$, $v_{n2}^{+}:={\stackrel{\dot{}}{y}}^{+}$, (103)

$v_{t1}^{-}:={\stackrel{\dot{}}{x}}^{-}+{f^{\prime }}_B{\stackrel{\dot{}}{y}}^{-}+r{\stackrel{\dot{}}{\phi }}^{-}$, $v_{t1}^{+}:={\stackrel{\dot{}}{x}}^{+}+{f^{\prime }}_B{\stackrel{\dot{}}{y}}^{+}+r{\stackrel{\dot{}}{\phi }}^{+}$, (104)

$v_{t2}^{-}:={\stackrel{\dot{}}{x}}^{-}+r{\stackrel{\dot{}}{\phi }}^{-}$, $v_{t2}^{+}:={\stackrel{\dot{}}{x}}^{+}+r{\stackrel{\dot{}}{\phi }}^{+}$. (105)

From the above-given equations we should determine the velocities of the wheels after the collision ${\stackrel{\dot{}}{x}}^{+}$, ${\stackrel{\dot{}}{y}}^{+}$, ${\stackrel{\dot{}}{\phi }}^{+}$ and reaction impulse multipliers ${\tilde{\lambda }}_1$, ${\tilde{\lambda }}_2$, ${\tau }_1$, ${\tau }_2$. Since it is a relatively complex system of nonlinear equations, an iterative method has been used to solve it. This method involves looking for a vector $Y:={\left[{\stackrel{\dot{}}{x}}^{+},{\stackrel{\dot{}}{y}}^{+},{\stackrel{\dot{}}{\phi }}^{+},{\tilde{\lambda }}_1,{\tilde{\lambda }}_2,{\tau }_1,{\tau }_2\right]}^{\text{T}}\in R^7$, which satisfies the equation

$Y=\Psi \left(Y\right)$, (106)

where $\Psi :R^7\to R^7$ denotes a mapping specified by the right sides of the above-given equations.

When we consider an elastic or elasto-plastic collision, then based on Equation (35) we should formulate relevant relations taking into account detailed data specifying the state of the wheel at point B. On this ground we establish the following equations:

${\tilde{\lambda }}_i^{\ast }=\beta {\tilde{\lambda }}_i,i=1,2$, (107)

$m{\stackrel{\dot{}}{x}}^{\ast }=m{\stackrel{\dot{}}{x}}^{+}-{f^{\prime }}_B{\tilde{\lambda }}_1^{\ast }-\mu \Vert G_1\Vert {\tilde{\lambda }}_1^{\ast }{\tau }_1^{\ast }-\mu {\tilde{\lambda }}_2^{\ast }{\tau }_2^{\ast }$, (108)

$m{\stackrel{\dot{}}{y}}^{\ast }=m{\stackrel{\dot{}}{y}}^{+}+{\tilde{\lambda }}_1^{\ast }-\mu {f^{\prime }}_B\Vert G_1\Vert {\tilde{\lambda }}_1^{\ast }{\tau }_1^{\ast }+{\tilde{\lambda }}_2^{\ast }$, (109)

$J{\stackrel{\dot{}}{\phi }}^{\ast }=J{\stackrel{\dot{}}{\phi }}^{+}-\mu \Vert G_1\Vert r{\tilde{\lambda }}_1^{\ast }{\tau }_1^{\ast }-\mu {\tilde{\lambda }}_2^{\ast }{\tau }_2^{\ast }$, (110)

${\tau }_i^{\ast }=\Pi \left({\tau }_i^{\ast }+\rho v_{ti}^{\ast }\right),\rho >0,i=1,2$, (111)

$v_{t1}^{\ast }:={\stackrel{\dot{}}{x}}^{\ast }+{f^{\prime }}_B{\stackrel{\dot{}}{y}}^{+}+r{\stackrel{\dot{}}{\phi }}^{\ast }$, (112)

$v_{t2}^{\ast }:={\stackrel{\dot{}}{x}}^{\ast }+r{\stackrel{\dot{}}{\phi }}^{\ast }$. (113)

On the basis of the above relations we can obtain two nonlinear equations in regard to ${\tau }_1^{\ast }$ and ${\tau }_2^{\ast }$, which we solve by iteration. Having determined their values we can determine the values of acceleration after the second phase of the considered collision ${\stackrel{\dot{}}{x}}^{\ast }$, ${\stackrel{\dot{}}{y}}^{\ast }$, ${\stackrel{\dot{}}{\phi }}^{\ast }$.

The relations given in this subsection refer to determining the acceleration of the wheel and its velocity after the collision when the wheel is positioned in segments IV, II and III and at point B (Figure 7). We do not present equations for segment I, as they are analogical to the equations for segment IV; while the equations at the singular point A, are similar to the above-described equations for point B.

6.3. Results of Computer Calculations

We consider planar motion of the non-deformable wheel on which constraints presented in Figure 1, Figure 2 and Figure 6 are acting. The wheel motion simulation is performed assuming the following values of the wheel’s parameters: the mass $m=10\text{kg}$, moment of inertia $J=0.05\text{kg}\cdot {\text{m}}^2$, and radius $r=0.1\text{m}$. The dimensions represented in Figure 1 are $h=5\text{cm}$ and $l=0.5\text{m}$. Furthermore, we assume that the gravitational acceleration is $g=9.81\text{m}/{\text{s}}^2$, and the initial conditions are defined with two vectors,

$X_{\text{o}}=\left[\begin{array}{c}0.8\text{m}\\ 0.2\text{m}\\ 0\text{rad}\end{array}\right]$, $V_{\text{o}}=\left[\begin{array}{c}-5\text{m}/\text{s}\\ -3\text{m}/\text{s}\\ 0\text{rad}/\text{s}\end{array}\right]$ (114)

Based on the relations presented in section 4, sets of relations specifying the acceleration of the wheel and forces of reaction and relations specifying discontinuous changes in the wheel’s velocity induced by the collision have been established for each segment of the boundary of constraints (Figure 6).

The simulation of the motion of the wheel has been performed over the time interval $t\in \left(0;1\right]\text{s}$. We assumed the value of the coefficient of friction between the wheel and the boundary of constraints, $\mu =0.4$. Two tasks differing only with the value of the coefficient of restitution $\beta$ have been considered. In the first task, the motion of the wheel has been simulated with the assumption of the hypothesis of plastic collisions $\beta =0$. In the case of the second task, we have assumed $\beta =0.3$, which corresponds to the elasto-plastic hypothesis.

The wheel position sequences during the simulation are illustrated in Figure 10 (plastic collisions) and Figure 11 (elasto-plastic collisions). The figures show also the initial conditions of the wheel. It follows from them that in the initial stage of its motion the wheel moves in a translational motion and after the first collision it starts rotating.

The next Figure 12 represents two trajectories (paths of the wheel in motion) followed by the centre point of the wheel in the two tasks under analysis.

The two last graphs illustrate the positions of the wheel’s centre point and the angle of rotation in time (Figure 13), as well as the speed of change of these values

Figure 10. The sequence of the positions of the wheel during simulation (plastic collisions).

Figure 11. The sequence of the positions of the wheel during simulation (elasto-plastic collisions).

Figure 12. Trajectories followed by the centre point of the wheel: plastic collisions (solid line), elasto-plastic collisions (dashed line).

Figure 13. Curves of wheel centre point positions (a) and (b) and of the rotation angle (c) over time: plastic collisions (solid line), elasto-plastic collisions (dashed line).

in time (Figure 14). The results illustrated in these figures reveal major differences in the behaviour of the wheel, depending on the collision hypothesis assumed.

7. Conclusions

An analysis of the non-smooth problem of the motion of a mechanical system in the form of a wheel with constrained motion has been presented in the paper. It has been demonstrated that the formulation of the non-smooth problem of motion encompasses relations determining accelerations of bodies and relations determining discontinuous changes in the velocity of the bodies induced by the

Figure 14. Curves of wheel centre point positions (a) and (b) and of the angular velocity (c) over time: plastic collisions (solid line), elasto-plastic collisions (dashed line).

collision. Descriptions of these relations have been established according to the hypotheses specifying interactions (reactions) between non-deformable bodies coming into contact.

In the wheel example under consideration, the interpretation of the assumed hypotheses is clear, and therefore the formulation of a non-smooth problem of motion is relatively simple. In the case of a rigid body in a three-dimensional space, the description of the reactions is more complex, and apart from that, it is impossible to unequivocally determine the acceleration and the reaction force for individual positions of the body. Such singular states of the system result from the approximate character of the hypotheses specifying the forces acting between the bodies in contact, which applies in particular to friction forces.

Non-smooth problems of motion expand on the body of classical problems of the theoretical mechanics, thus increasing the scope of modeling methods and simulation studies on the motion of and loads applied to machines, vehicles and components of engineering structures [15].

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Grzesikiewicz, W. and Zbiciak, A. (2018) Opis ruchu układu mechanicznego z więzami jednostronnymi (Description of Motion of Mechanical System with Unilateral Constraints). Oficyna Wydawnicza Politechniki Warszawskiej, Warszawa.
[2]	Painlevé, P. (1895) Leçons sur le frottement. Hermann, Paris.
[3]	Routh, E.J. (1905) Dynamics of a System of Rigid Bodies. MacMillan, London.
[4]	Jean, M. and Moreau, J.J. (1992) Unilaterality and Dry Friction in the Dynamics of Rigid Body Collections. Proceedings of Contact Mechanics International Symposium, Lausanne, 31-48.
[5]	Jean, M. (1999) The Nonsmooth Contact Dynamics Method. Computer Methods in Applied Mechanics and Engineering, 177, 235-257. https://doi.org/10.1016/S0045-7825(98)00383-1
[6]	Moreau, J.J. and Panagiotopoulos, P.D. (1988) Non-Smooth Mechanics and Applications. In: CISM Courses and Lectures, Vol. 302, Springer, Wien.
[7]	Panagiotopoulos, P.D. (1985) Inequality Problems in Mechanics and Applications: Convex and Nonconvex Energy Functions. Birkhäuser, Basel.
[8]	Panagiotopoulos, P.D. (1993) Hemivariational Inequalities: Applications in Mechanics and Engineering. Springer-Verlag, Berlin. https://doi.org/10.1007/978-3-642-51677-1
[9]	Grzesikiewicz, W. (1990) Dynamika układów mechanicznych z więzami (Dynamics of Mechanical Systems with Constraints). Prace Naukowe Politechniki Warszawskiej, Mechanika z. 117, Wydawnictwa Politechniki Warszawskiej, Warszawa.
[10]	Grzesikiewicz, W. and Wakulicz, A. (1999) Dynamics of Non-Smooth Mechanical System. Machine Dynamics Problems, 23, 25-37.
[11]	Grzesikiewicz, W., Wakulicz, A. and Zbiciak, A. (2009) Succession of Constraint Imposed on Time Function. In: Monograph: Polioptimization and Computer Aided Design, Vol. 7, Publishing House of Koszalin University of Technology, Koszalin, 49-56.
[12]	Stronge, W.J. (2013) Smooth Dynamics of Oblique Impact with Friction. International Journal of Impact Engineering, 51, 36-49. https://doi.org/10.1016/j.ijimpeng.2012.08.001
[13]	Stewart, D.E. (2000) Rigid-Body Dynamics with Friction and Impact. SIAM Review, Society for Industrial and Applied Mathematics, 42, 3-39. https://doi.org/10.1137/S0036144599360110
[14]	Le, X.A. (2003) Dynamics of Mechanical Systems with Coulomb Friction. Springer-Verlag, Berlin, Heidelberg.
[15]	Dubois, F., Acary, V. and Jean, M. (2018) The Contact Dynamics Method: A Nonsmooth Story. Comptes Rendus Mécanique, 346, 247-262. https://doi.org/10.1016/j.crme.2017.12.009