Non-smooth dynamic modeling and simulation of an unmanned bicycle on a curved pavement

The non-smooth dynamic model of an unmanned bicycle is established to study the contact-separate and stick-slip non-smooth phenomena between wheels and the ground. According to the Carvallo-Whipple configuration, the unmanned bicycle is reduced to four rigid bodies, namely, rear wheel, rear frame, front fork, and front wheel, which are connected by perfect revolute joints. The interaction between each wheel and the ground is simplified as the normal contact force and the friction force at the contact point, and these forces are described by the Hunt-Crossley contact force model and the LuGre friction force model, respectively. According to the characteristics of flat and curved pavements, calculation methods for contact forces and their generalized forces are presented. The dynamics of the system is modeled by the Lagrange equations of the first kind, a numerical solution algorithm of the dynamic equations is presented, and the Baumgarte stabilization method is used to restrict the drift of the constraints. The correctness of the dynamic model and the numerical algorithm is verified in comparison with the previous studies. The feasibility of the proposed model is demonstrated by simulations under different motion states.


Introduction
Bicycles, after more than 200 years of their development, have become the most widely used transportation tools in the world. Despite their simple mechanical structure, bicycles exhibit complex dynamical behaviors. Although bicycles and motorcycles have certain differences, their dynamics and control problems are very similar, and thus we do not explicitly distinguish between them. Researchers have extensively investigated the dynamic modeling [1][2][3] , stability analysis [4][5] , path tracking [6][7][8] , and equilibrium control [9][10][11] of bicycles. Carvallo [12] and Whipple [13] first simplified bicycles into multibody systems consisting of four rigid bodies, developed a model of bicycle dynamics, and linearized the resulting dynamic equations. Getz and Marsden [14] developed a simplified dynamic model of bicycles according to their geometry, mass, and constraints. Meijaard et al. [15] extended the Carvallo-Whipple bicycle model by considering tire geometry, road slope, the damping between tires and the ground, and air resistance. Xiong et al. [16][17] studied the dynamics and stability of an idealized benchmark bicycle moving on a revolute surface. Basu-Mandal et al. [18] developed a bicycle dynamic model based on the Lagrange equations and the Newton-Euler method and analyzed the circular motion of an unmanned bicycle. Pappalardo et al. [19][20] investigated the linear and nonlinear stability of multi-body systems with complete and incomplete constraints and discussed the effects of different parameters on bicycle stability. Chen et al. [21][22] established a 9-DOF (degrees of freedom) bicycle dynamic model based on the Lagrange equations and calculated the equilibrium point of bicycle movement. The self-stability of bicycles was studied and designed a fuzzy controller to keep bicycles in balance [23][24] . Keo and Yamakita [10] established a simplified bicycle model by ignoring the mass and inertia of wheels and installed balancers on bicycles to facilitate balance control.
Extensive research has been conducted on the dynamic modeling and control of bicycles. However, these studies assume that wheels experience pure rolling on the ground and describe the interactions between wheels and the ground as constraints, such as permanent contact constraints and pure rolling constraints, without considering the slipping of wheels relative to the ground and their separation from the ground. These assumptions simplify the dynamic modeling of bicycles and do not reflect the non-smooth events between wheels and the ground, such as contact-impact [25] and stick-slip [26] . Unmanned bicycles cannot avoid non-smooth events while doing head-lifting, flying, drifting, and other large maneuvers, and thus non-smooth dynamic modeling is the basis for realizing large maneuvers of unmanned bicycles. Sharp [27] formulated linearized equations for bicycles considering the wheel slipping effect. Hauser and Saccon [28] studied the extreme motion state of motorcycles in which wheels do not contact the ground. Yi et al. [29][30] established a dynamic model for an unmanned motorcycle by considering the slipping between wheels and the ground. Gerdes's group [31][32] established a dynamic model for the drift maneuver of unmanned vehicles and conducted a series of experiments on drift control. However, non-smooth events of contact-separate and stick-slip between the bicycle wheels and the ground are not comprehensively studied in these analyses, and driving problems on curved pavements are rarely investigated [33] .
Over the past two decades, multifarious numerical algorithms and experimental studies have been presented to analyze the non-smooth dynamics of multibody systems [34][35][36][37][38][39] . Non-smooth contact forces mainly include normal and tangential (friction) forces. To study the non-smooth dynamics of unmanned bicycles, it is necessary to establish a suitable mechanical model for contact-impact-friction problems. The Hunt-Crossley contact force model proposed by Hunt and Crossley [40] is the most widely used normal contact force model, and it can describe contact problems more rationally than the Hertz contact force model [41] and the linear spring-damping model [42] because of its ability to describe nonlinear characteristics and energy loss during contact. The Coulomb dry friction model is the most widely known friction model, and it can describe stick-slip motion effectively. However, when the relative velocity of contact points is zero, the friction force is a multi-valued function of acceleration, causing great difficulties for numerical calculations [26,43] . The modified Coulomb friction model expresses friction as a single-valued function of velocity; however, it fails to describe stick-slip motion [37] . The LuGre friction model proposed by Canudas et al. [44] can successfully describe stick-slip motion and the Stribeck effect; however, the difficulty of obtaining physical parameters limits the use of this model. In our previous research [45] , a parameter identification method for the LuGre friction model was proposed and then tested on a wheel mechanism. There were also many tire modeling and tire force models proposed [46][47][48][49][50][51] ; however, their main concern is the modeling of tire sideslip, rather than the modeling of stick-slip friction.
In the present work, the non-smooth dynamic model of an unmanned bicycle moving on a curved pavement is established based on the Hunt-Crossley contact force model and the LuGre friction model to study the contact-separate and stick-slip non-smooth phenomena between wheels and the ground. The remaining article is organized as follows. The system description of the unmanned bicycle and the adopted contact force models are presented in Section 2. The generalized forces of the wheel-ground contact forces for the unmanned bicycle are derived in Section 3. The dynamic model of the unmanned bicycle and its solution method are presented in Section 4. The effectiveness of the dynamic modeling method and the numerical algorithm is verified by several numerical examples in Section 5. In Section 6, the main conclusions are drawn.
2 Unmanned bicycle system and contact force models 2.1 Unmanned bicycle system The Carvallo-Whipple configuration shown in Fig. 1 is the most commonly used bicycle configuration. The bicycle consists of four rigid bodies, namely, rear wheel, rear frame, front fork, and front wheel (numbered as 1, 2, 3, 4, respectively). For convenience, the inertial reference frame Oxyz is established on the ground, and the body-fixed reference frame C i x i y i z i (i = 1, 2, 3, 4) is defined on the rigid body i (C i is the center of mass of the corresponding rigid body). Suppose that the mass and inertia matrix of each rigid body with respect to the body-fixed reference frame are m i and J i , respectively. In addition, the following assumptions are made in this analysis.
(I) The rigid bodies of the bicycle are connected by perfect revolute joints, and the clearance and friction in the joints are ignored.
(II) The front and rear wheels have a knife-edge construction, and the contact between the wheels and the ground is regarded as point contact.
(III) The front and rear wheels are circumferentially symmetrical (the center of mass coincides with the center of the circle).
(IV) The rear frame and the front fork are symmetrical with respect to the plane of the rear and front wheels.
Geometric parameters are assigned when the unmanned bicycle is upright on the ground and the steering angle, i.e., the rotation angle of the front fork relative to the rear frame is zero. In Fig. 2, the radii of the rear and front wheels are R 1 and R 4 , respectively, and the positions of the centers of mass of the rear frame and the front fork are (y 2 , z 2 ) and (y 3 , z 3 ), respectively. Let the distance between the contact points of the front and rear wheels be w, the distance of the front wheel-ground contact point after the intersection of the steering axis and the ground (the trail) be c, and the angle between the steering axis and the vertical line be η.
In this analysis, the mass center coordinates and the Euler angles of the rigid bodies are used to describe their positions and orientations; hence, the generalized coordinates of the unmanned bicycle can be expressed as

Rear wheel
Rear frame Steering Front fork Front wheel Fig. 1 Schematic representation of the unmanned bicycle (color online) where in which (x i , y i , z i ) are the position coordinates of the mass center C i of the rigid body i in the inertial reference frame, and (ψ i , θ i , ϕ i ) are the Euler angles of the body-fixed frame C i x i y i z i with respect to the inertial frame. Without loss of generality, the Euler angle rotation order of 3-1-3 is used in this analysis, as shown in Fig. 3. Here, we define that the tilt and pitch angles of the rear frame are θ 2 − π/2 and ϕ 2 , respectively.

Normal contact force model
When the unmanned bicycle wheels contact the ground, normal and tangential (friction) contact forces will be generated; thus, suitable contact force models should be selected to calculate these contact forces. In the Hunt-Crossley contact force model [40] , when two objects are in contact, the normal contact force is expressed as a function of the penetration depth and the penetration velocity, and this model is widely used since it can accurately describe the damping and nonlinear characteristics of a contact process. The Hunt-Crossley contact force model can be described as where δ i andδ i represent the penetration depth and penetration velocity, respectively, between wheel i (i = 1, 4) and the ground. K is the stiffness coefficient, χ is the hysteresis factor, and the values of m and n are related to the material and geometrical properties of the contact area, respectively. The Hunt-Crossley contact force model is illustrated in Fig. 4. Numerous methods have been proposed to determine the hysteresis factor of the Hunt-Crossley contact force model [52] . However, as no unified conclusion is yet available, the classical Hunt-Crossley contact force model is used in this analysis.

Tangential contact force model
The Coulomb dry friction model is generally used to calculate friction between two objects. This model can well describe static friction characteristics and the stick-slip motion between two objects; however, it is difficult to judge the stick-slip state when the relative tangential velocity of the contact points is zero. The LuGre friction model is the most comprehensive model to calculate friction between two objects, and it can describe a variety of phenomena including stick-slip motion, the Stribeck effect, viscous friction, pre-slipping, and friction hysteresis. This model also has continuity, which creates convenience to numerical calculations. Therefore, the LuGre friction model is used in this analysis to calculate tangential contact forces between the wheels and the pavement. The LuGre friction model can be described as follows [44] : where z is the average deformation of bristles, as shown in Fig. 5,ż is the relative velocity between two surfaces, σ 0 is the bristle stiffness coefficient, σ 1 is the micro-damping coefficient, σ 2 is the viscous coefficient, v r is the vector of relative tangential velocity, v s is the Stribeck velocity, g(v r ) is a velocity-dependent function that can reproduce both the Coulomb friction and the Stribeck effect, µ c is the kinetic friction coefficient, and µ s is the static friction coefficient. The main disadvantage of the LuGre friction model is the difficulty in the selection of the parameters σ 0 and σ 1 . In our previous study [45] , a suitable method for the identification of these two parameters was proposed.
3 Generalized forces of contact forces between wheels and the ground 3.1 Flat pavement driving problems When the unmanned bicycle moves on a flat pavement, the calculation of contact forces is relatively simple, as the normal contact force is vertical and the friction force is located in the horizontal plane. For the rigid body i (i = 1, 4), the penetration depth between the wheel and the ground can be expressed by the generalized coordinates of the system, i.e., Taking the first derivative of Eq. (5) with respect to time yieldṡ The normal contact force between the wheels and the pavement can be calculated by substituting Eqs. (5) and (6) into Eq. (3). Let the coordinates of the wheel-ground contact point P i in the inertial reference frame be P i (x Pi , y Pi , z Pi ); hence, the generalized force vector of the normal contact force can be expressed as where Therefore, the generalized force vector of the normal contact force between the wheels and the pavement is As shown in Fig. 6, the relative tangential velocity of wheel-ground contact points can be divided into velocity v τ i along the wheel rolling direction and velocity v ni perpendicular to the wheel rolling direction, i.e., Fig. 6 Relative velocity of contact points (color online) Therefore, the relative tangential velocity of wheel-ground contact points can be expressed as where v xi and v yi are the components of tangential velocity v ri along the x-and y-axes of the inertial frame, respectively. Now, the friction force F fi between the wheels and the pavement can be calculated by substituting Eq. (11) into Eq. (4). The generalized force vector of friction forces can be expressed as where Therefore, the generalized force vector of friction forces between the wheels and the pavement is 3.2 Driving scenario on a curved pavement 3.2.1 Definition of a curved pavement The driving of an unmanned bicycle on a curved pavement surface is relatively complicated. The curved road surface shown in Fig. 7 can be defined as Let e 1 , e 2 , and e 3 be unit vectors of the inertial reference frame. For any potential wheelground contact point P (x, y, f (x, y)) on the surface, denote Thus, a set of orthogonal unit vectors is obtained, where g 1 and g 2 determine the tangential planes of the curved surface at the contact point, and g 3 is the normal vector of the plane at the contact point.

Determination of wheel-ground contact points
The Hunt-Crossley contact force model assumes a certain amount of penetration between wheels and the pavement, generating a contact area rather than a contact point. To calculate wheel-ground contact forces, it is first necessary to determine wheel-ground contact points, which are considered as points with the largest penetration depth in the contact area. In Fig. 8, the point closest to the wheel center on the intersection line between the wheel i (i = 1, 4) and the curved pavement is set as P i . The extension line connecting C i and P i intersects the wheel edge at P i . Thus, the contact points on the pavement and the wheel are P i and P i , respectively. Let the coordinates of the wheel center C i be (x Ci , y Ci , z Ci ), and the plane defined by the wheel be S i ; hence, the constrained extremum problem can be used to calculate the position coordinates (x Pi , y Pi , z Pi ) of point P i , which is described as Equation (17) can be solved by the fmincon function in MATLAB. The coordinates of point P i can be determined by C i and P i . Suppose that the relative position of point P i on the wheel is denoted as γ i , and thus the vector diameter r CiP i in the body-fixed reference frame of the wheel can be expressed as r CiP i = (cos γ i , sin γ i , 0) T .

Generalized forces of contact forces
The wheel-ground penetration depth δ i and the penetration velocityδ i can be expressed as where v P i is the velocity of point P i . Let the velocity of the wheel center C i be v Ci , and the angular velocity of the wheel i in the body-fixed reference frame be ω i . Therefore, in which D i is the rotation matrix of rigid body i. The wheel-ground normal contact force F ni along the direction of g 3 can be calculated by substituting Eq. (18) into Eq. (3), which yields The generalized force vector can be expressed as where Therefore, the generalized force vector of the normal contact force between the wheels and the pavement can also be expressed in the form of Eq. (9).
The relative tangential velocity of wheel-ground contact points can be obtained by projecting velocity v P i to the tangent vectors g 1 and g 2 , i.e., The friction force on the wheels can be obtained by substituting Eq. (23) into Eq. (4), which yields The generalized force vector of friction force can be determined as where W i is the same as that in Eq. (22). Therefore, the generalized force vector of friction force between the wheels and the pavement is also expressed as Eq. (14). The limitation of the above wheel-ground contact force calculation process is only noticed when the first-order partial derivatives exist in the surface equation of the pavement, and each wheel has only one contact point with the pavement; however, many actual situations can be considered to satisfy its applicable conditions.

Derivation and solution of dynamic equations of the unmanned bicycle 4.1 Constraint equations of revolute joints
As shown in Fig. 9, the rigid bodies of the unmanned bicycle system are connected by perfect revolute joints, which restrict the rotation of the rigid bodies only around their axis. These joints require that point A on the rotating shaft be the same point on rigid bodies i and j, namely, where r Ci and r Cj are the position vectors pointing to points C i and C j from the origin of the inertial frame, respectively, r CiA and r Cj A are the vectors in the corresponding body-fixed reference frame expressions pointing to point A from points C i and C j , respectively, and D i and D j are the rotation matrices, similar to that of Eq. (19), of rigid bodies i and j, respectively. In addition, for these revolute joints, the rotational axes on rigid bodies i and j must coincide. Hence, where l i ij and l j ij are the direction cosine vectors of the rotation axis in the corresponding bodyfixed frames. Only two of the three scalars of Eq. (27) are independent, and they are denoted by Φ ij2 ; hence, the revolute joint between rigid bodies i and j produces five constraints as Therefore, the constraints of the unmanned bicycle system can be expressed as The above constraints are all holonomic. There are 3 × 5 = 15 constraint equations in total in the revolute joints, and thus the DOF of the unmanned bicycle 4 × 6 − 15 = 9.

Derivation and solution of dynamic equations
According to the Lagrange equations of the first kind, the dynamic equations of the unmanned bicycle can be expressed as where T is the kinetic energy of the system, Q is the generalized force vector of non-contact forces, including gravity, rear wheel driving moment, and steering moment, Q n and Q f are generalized force vectors of contact forces given by Eqs. (9) and (14), respectively, Φ is the constraint vector of revolute joints shown in Eq. (29), Φ q = ∂Φ ∂q , and λ is the corresponding Lagrange multiplier vector.

The kinetic energy of the unmanned bicycle is
where M is the mass matrix of the bicycle system, in which The Baumgarte constraint stabilization method [53] is used to restrict the drift of the constraint equationsΦ The expansion of the above equation leads to whereΦ Now, substituting Eq. (31) into Eq. (30) results in The vector of the Lagrange multipliers can be obtained by substituting Eq. (37) into Eq. (35), which yields The differential equations of motion for the unmanned bicycle can be obtained by substituting Eq. (39) into Eq. (37), which results in The above equation can be solved by using any common numerical algorithm for ordinary differential equations [54] .

Numerical examples
Meijaard et al. [15] defined a benchmark bicycle. To verify the effectiveness of the dynamic modeling method and the numerical algorithm proposed in this analysis and to investigate non-smooth dynamics (contact-separate and stick-slip) between the unmanned bicycle wheels and the pavement, the mass and geometric parameters used for the benchmark bicycle are selected for simulations, and simulation results are compared with the findings of Schwab and Meijaard [55] . The parameters of the benchmark bicycle are presented in Table 1.

Self-stability characteristics of the unmanned bicycle
Self-stability is the ability of an unmanned bicycle to automatically recover from lateral disturbances and maintain its balance under unmanned operation. Self-stability is one of the important characteristics of an unmanned bicycle. Schwab and Meijaard [55] pointed out that when the bicycle velocity was low, the swing amplitude increased gradually after lateral disturbances, and when the driving velocity was high, the bicycle experienced a spiral motion with a gradually decreasing radius of curvature after lateral disturbances. It was found that the bicycles recovered from lateral interferences with stability within a certain velocity range. When the forward velocity (v) of the benchmark bicycle was 4.29 m/s < v < 6.02 m/s, it was asymptotically stable. In this paper, the self-stability of the unmanned bicycles is examined under three different initial forward velocities of 3.8 m/s, 4.5 m/s, and 6.8 m/s. It should be noted that the initial tilt angle, the steering angle, and the corresponding angular velocities of the unmanned bicycle are all zero.
The initial velocity of the unmanned bicycle is set as 4.5 m/s, which is in the asymptotic stability interval proposed by Schwab and Meijaard [55] . When t = 2 s, an interference force of 500 N is applied at the mass center of the rear frame perpendicularly to the frame plane for 0.1 s. Figures 10(a) and 10(b) display the time histories of the tilt and steering angles of our multibody model compared with Meijaard's linear model [55] . The phase diagrams of the tilt and steering angles of the bicycle multibody model are also presented in Figs. 10(c) and 10(d). It is noticeable that the tilt and steering angles quickly deviate from zero after the bicycle is disturbed and gradually converge to zero after certain fluctuations, reflecting that the bicycle has the characteristics of asymptotic stability, which is consistent with Meijaard's conclusion [55] . Meijaard's linear model [55] Present multibody model Meijaard's linear model [55]  The results of the Baumgrate stabilization method are exhibited in Fig. 11. It is noticeable that when the Baumgrate stabilization method (ξ B = 100 and ζ B = 100) is used to restrict the constraint drift, the constraint drift only increases to a magnitude of 10 −8 when the interference suddenly occurs and then rapidly decreases to a magnitude of 10 −10 . Without the Baumgrate stabilization method (ξ B = 0 and ζ B = 0), the constraint drift increases rapidly to a magnitude of 10 −5 and then manifests a growing trend, indicating that the Baumgrate stabilization method has effective influence on controlling the constraint drift problems.
The initial velocity of the unmanned bicycle is set as 3.8 m/s. When t = 2 s, a lateral disturbance force of 0.5 N perpendicular to the rear frame plane is applied at the mass center for 0.1 s. Figure 12 displays the time histories of the tilt angle and the steering angle of two bicycle models. The phase diagrams of the tilt and steering angles of our multibody model are also presented in Fig. 12. It is noticeable from Figs. 12(a) and 12(b) that when the bicycle is slightly disturbed, both its tilt and steering angles start to fluctuate, and the amplitude increases rapidly. When the tilt and steering angles increase to a certain degree, the role of nonlinear terms in the dynamic equations cannot be neglected, and they finally enter a state of periodic variation, which is expressed by limit loops in the phase diagrams shown in Figs. 12(c) and 12(d). However, the linear model proposed in Ref. [55] is no longer applicable at this time.
The initial velocity of the unmanned bicycle is set as 6.8 m/s. When t = 2 s, a lateral disturbance force of 500 N is applied for 0.1 s. The time histories of the tilt and the steering angles of our multibody model compared with Meijaard's linear model are exhibited in Fig. 13.  Steering angle/rad (d)

Present multibody model
Meijaard's linear model [55] Present multibody model Meijaard's linear model [55]  Meijaard's linear model [55] Present multibody model Meijaard's linear model [55] Time/s (a) It is evident that after the unmanned bicycle is subject to lateral interferences, the tilt angle and the steering angle of the unmanned bicycle do not converge after a period of fluctuation and increase gradually, manifesting the characteristics of instability. The asymptotically stable, non-asymptotically stable, and unstable motions of the unmanned bicycle can be noticed in the above three cases, and these results are consistent with the analytical results in Ref. [55], verifying the correctness of the proposed multibody model of the unmanned bicycle.

Rear wheel slipping movement of the unmanned bicycle
The rear wheel slipping or head lifting movement of the unmanned bicycle can be achieved by applying a suitable driving moment on the rear wheel. In this part, the rear wheel driving moment is set as M 1 = 100e −(t−5) 2 N · m, and the initial forward velocity of the unmanned bicycle is 5.4 m/s. Figure 14 presents the time histories of the relative tangential velocities, normal contact forces, and friction forces at wheel-ground contact points, and the angular velocities of front and rear wheels. It is noticeable that, with the increase of the driving moment, the unmanned bicycle accelerates and the rear wheel friction force gradually increases. Due to the inertia force and the reaction driving moment from the rear wheel to the rear frame, the wheel-ground normal contact force of the rear wheel gradually increases, whereas the front wheel contact force gradually decreases. When the driving moment increases to a certain degree, the friction force is not sufficient to support the pure rolling of the rear wheel; hence, slipping occurs between the rear wheel and the pavement. At this moment, the rear wheel friction changes from static friction to kinetic friction. As the driving moment decreases, the relative tangential velocity of rear wheel-ground contact points gradually decreases to zero again, and the friction changes back to static friction. It should be noted that the front wheel always maintains a pure rolling state and never slips.

Head lifting movement of the unmanned bicycle
In this part, a head lifting movement is achieved by applying an appropriate driving moment on the rear wheel of the unmanned bicycle. The driving moment is set as M 1 = 134e −15(t−2) 2 N·m, and the initial forward velocity of the unmanned bicycle is 3 m/s. Figure 15 displays the time histories of the rear frame pitch angle, the angular velocities of front and rear wheels, and the normal contact forces and friction forces at wheel-ground contact points. With the increase in the driving moment, the unmanned bicycle accelerates. When the acceleration is large enough, owing to the combined effect of the inertia force and reaction driving moment from the rear wheel to the rear frame, the unmanned bicycle lifts its head. At this moment, the front wheel leaves the ground, and the contact force of the front wheel on the ground becomes zero. As the driving moment decreases, the front wheel again touches the ground. At the moment of impact between the front wheel and the ground, the normal contact force and the friction force change rapidly and then quickly become stable.

Downward motion of the unmanned bicycle on a stepped pavement
In this part, the problem of the unmanned bicycle driving on a stepped pavement is examined. As shown in Fig. 16, the initial forward velocity of the unmanned bicycle is 5 m/s, and the unmanned bicycle jumps down the steps twice. The time histories of the rear frame pitch angle, the angular velocities of the front and rear wheels, and the normal contact forces between the wheels and the ground are displayed in Fig. 17. It is noticeable that a contact-impact phenomenon occurs between the wheels and the pavement. In both drops, the front wheel of the bicycle falls first followed by the rear wheel, and there is a short time when both wheels leave the ground simultaneously. At the moment of collision between the wheels and the ground, a large impact force is generated, causing a large instantaneous peak of the friction force and then a sudden change in the angular velocities of the wheels. ) Pitch angle/rad x , x 0, where the unit of z is m. The initial velocity of the unmanned bicycle is 5 m/s. Figure 19 displays the time histories of the normal contact forces and the friction forces at wheel-ground contact points. It can be seen that the normal contact forces change suddenly, when the unmanned bicycle is driven from the flat surface to the cosine-shaped curved pavement. Both the normal contact force and the friction force change continuously with the terrain. The normal contact forces change regularly, as the largest contact forces are measured at the lowest points of the pavement surface. Therefore, it can be inferred that the proposed dynamic model can effectively simulate the driving problem of the unmanned bicycle on curved pavement surfaces.

Conclusions
In this paper, the dynamic model of an unmanned bicycle is established, and a numerical solution algorithm is presented based on the non-smooth dynamics of a multibody system. The contact between wheels and the ground is regarded as point contact, and the interaction between each wheel and the ground is simplified as the normal contact force and the friction force at the contact point. The Hunt-Crossley contact force model is used to describe the normal contact force, and the LuGre friction model is adopted to calculate the friction force. The nonsmooth dynamic model of the bicycle is developed based on the Lagrange equations of the first kind, and the Baumgrate constraint stabilization method is used to effectively restrict the drift of constraint equations during numerical calculations. The determination of the contact points between the wheels and the curved pavement is transformed into a constrained extremum problem. The effectiveness of the multibody dynamic model is compared with the research work of Meijaard. Finally, several numerical examples are presented to simulate the bicycle rear wheel slipping motion, the bicycle head lifting motion, the driving scenario on a stepped pavement, and the driving scenario on a curved pavement. It is demonstrated that the proposed dynamic model and the numerical algorithm effectively reflect the non-smooth dynamic behaviors of the unmanned bicycle.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.