A Multibody System Approach for the Systematic Development of a Closed-Chain Kinematic Model for Two-Wheeled Vehicles

Abstract

In this investigation, a closed-chain kinematic model for two-wheeled vehicles is devised. The kinematic model developed in this work is general and, therefore, it is suitable for describing the complex geometry of the motion of both bicycles and motorcycles. Since the proposed kinematic model is systematically developed in the paper by employing a sound multibody system approach, which is grounded on the use of a straightforward closed-chain kinematic description, it allows for readily evaluating the effectiveness of two alternative methods to formulate the wheel-road contact constraints. The methods employed for this purpose are a technique based on the geometry of the vector cross-product and a strategy based on a simple surface parameterization of the front wheel. To this end, considering a kinematically driven vehicle system, a comparative analysis is performed to analyze the geometry of the contact between the front wheel of the vehicle and the ground, which represents a fundamental problem in the study of the motion of two-wheeled vehicles in general. Subsequently, an exhaustive and extensive numerical analysis, based on the systematic multibody approach mentioned before, is carried out in this work to study the system kinematics in detail. Furthermore, the orientation of the front assembly, which includes the frontal fork, the handlebars, and the front wheel in a seamless subsystem, is implicitly formulated through the definition of three successive rotations, and this approach is used to propose an explicit formulation of its inherent set of Euler angles. In general, the numerical results developed in the present work compare favorably with those found in the literature about vehicle kinematics and contact geometry.

1. Introduction

This paper investigates the kinematic modeling of two-wheeled vehicles, paying particular attention to the contact constraints of the wheels with the road profile. As a distinguishing feature of this study, which is based on a fully nonlinear multibody description of the system geometry, two methods for the formulation of the contact constraints are analyzed. This is done without recurring to the common linearization strategy found in the literature. Additionally, considering a closed-chain multibody approach for modeling the vehicle kinematics, the angular orientation of the front assembly (the frontal fork, the handlebars, and the front wheel) is implicitly and explicitly formulated in this work. In the remaining parts of this section, the background of the research topic, the definition of the problem of interest for this research work, a discussion of the literature concerning the issues at hand, the contributions of the present research, and the organization of the entire manuscript are reported.

1.1. Background Information and Research Significance

The mathematical modeling of two-wheeled vehicles is a well-established research topic dating back to the end of the 19th century when the first relevant mathematical model was proposed [1]. Since then, research about two-wheeled vehicle dynamics has been evolving. However, it is still challenging to rigorously describe the complex motion of a two-wheeled vehicle [2], and many questions still need to be resolved [3]. For instance, despite the apparent simplicity of the topological structure of this mechanical system, there is no complete model capable of exhibiting and capturing its more representative dynamical aspects in a simple and manageable manner. Thus, unless one considers advanced multibody models for detailed numerical simulations [4], which offer little insight into the physical phenomena due to their inherent complexity, the kinematic and dynamic modeling of two-wheeled vehicles is still an open issue. It is, therefore, desirable to develop in a parametric fashion a comprehensive mathematical model of a general two-wheeled vehicle to try to combine the benefits of the simplified physical-based linear models and the systematic nonlinear multibody approach. Reaching this goal will allow for capturing the fundamental kinematic and dynamic aspects of this family of mechanical systems by performing a trade-off between a manageable degree of complexity and high accuracy.

The initial interest of the scientific community in two-wheeled vehicles arose mainly from the search for an explanation for the self-stability experimentally observed in bicycles. This phenomenon consists of the fact that a bicycle is statically unstable. Still, it can automatically self-steer in a specific forward speed range, recovering from relatively small perturbations that would induce the fall, thereby exhibiting a certain type of dynamic stability [5]. For instance, the so-called Whipple–Carvallo bicycle model predicted a range of linear velocity for the bicycle self-stability confirmed by practical experiments [1,6]. At first, it was considered that this was due to the gyroscopic precession of the bicycle front wheel. However, a thorough discussion of Kooijman et al. provided a more comprehensive explanation, reporting the possibility of bicycle self-stability without this gyroscopic effect [7]. Subsequently, the interest of motorcycle producers in the mid-1950s stimulated the investigation of more complex and realistic physical phenomena found in two-wheeled vehicles in general. Fundamental research works by Sharp [8], among other researchers [9,10,11], to name a few, made it possible to identify the principal vibration modes of motorcycles and their leading causes for assessing vehicle stability. In particular, Sharp computed the principal eigenvalues of the linearized model of a motorcycle, parameterized again by the vehicle forward velocity [8]. In this way, the three main vibration modes of two-wheeled vehicles were identified and named as the wobble, the weave, and the capsize, respectively. Thus, Sharp’s research work is considered the starting point of modern modeling of motorcycle dynamics [12].

Due to the research emphasis on accurately reproducing realistic physical effects, mathematical models of two-wheeled vehicles were developed with as much complexity as manual calculation allowed. However, the numerical solution of the equations of motion was not an option at the time, so the linearization approach and subsequent stability analysis by the eigenvalues criterion was the central methodology used. The gradual appearance of computational tools would open the doors to the systematic development of complex models using redundant coordinate approaches typical of advanced multibody dynamics computational tools. This allowed the use of numerical solutions for kinematic and dynamic models with a high level of detail, such as those developed by Sharp [13] and Cossalter [14], to name a few. Further refinements led to the advanced mathematical models currently available in the literature [15,16], and some of them are based on the multibody approach to vehicle dynamics [17,18,19,20,21]. These are capable of simulating the behavior of two-wheeled vehicles with high fidelity, being able to reproduce instabilities under operating conditions based on data recorded during field operations on the actual vehicles. For instance, the reproduction through numerical simulation of the chattering behavior [22,23], the shimmy mode [4,12], or the wheel patter instability [24] showed a high degree of correlation between the numerical and the experimental results, thereby demonstrating the considerable level of maturity reached by the research on this topic.

1.2. Formulation of the Problem of Interest for This Investigation

The number of two-wheeled vehicles in use is constantly increasing. In particular, motorcycles are an essential mode of transportation, and this trend, combined with the growing popularity of eBikes, shows that these vehicles are part of future transportation [25,26]. Two-wheeled vehicles, however, have a high accident rate. The evidence indicates that they are among the most dangerous modes of transport [27,28,29,30]. Alternatives to address this problem include at least two interventions: (a) the adjustment and improvement of the handling qualities of two-wheeled vehicles; and (b) the inclusion of passive and active safety features in two-wheeled vehicles. In particular, the fact that human errors are one of the most contributing factors to the causation of accidents due to loss of control is highlighted in the literature [31]. Such a situation occurs in 32% of accidents involving two-wheeled vehicles [32]. Therefore, it is of interest to add two-wheeled vehicles the safety features implemented in cars, such as anti-lock braking, traction control, or stability systems, given that such systems resulted in an excellent accident reduction for cars [33,34,35,36]. Thus, adequate systems could help the rider stay in control of the two-wheeled vehicle during dangerous situations and avoid accidents [37,38]. However, the same results do not occur when applying such systems for two-wheeled vehicles [39,40]. The main reason for these inefficiencies lies in the incongruity in the two-wheeled vehicle models used to develop them, the understanding of their dynamics, the simulation of the driver inputs, and the impact of the infrastructure characteristics [27]. In general, the efficiency related to the model-based approach is often closely linked to their capability to be faithful regarding the actual system dynamics in conjunction with the intervention of the vehicle driver. Furthermore, vehicle instability is one limiting reason that automated emergency braking, collision avoidance systems, and similar active safety technologies have not been developed for motorcycles and two-wheeled vehicles in general [41,42].

In the design and analysis of two-wheeled vehicle systems, it is clear that a better understanding of vehicle dynamics variables will improve the effectiveness of driving assistance systems [43,44]. For this reason, many researchers have used the linearized benchmark model based on the Whipple–Carvallo bicycle system [45]. This benchmark model comprises four rigid bodies connected by revolute joints, including the wheels making a knife-edge rolling point contact with the road [46]. Despite its simplifying assumptions, the Whipple–Carvallo bicycle model can correctly reproduce the stability features of a wide family of two-wheeled vehicle systems in good agreement with the experimental results. However, linear models only validate a limited range of parameter values and ignore some features such as zero velocity balancing or hand-free curve riding [47,48]. Besides, to effectively deploy electronic control systems for two-wheeled vehicles, it is required to know the instantaneous dynamic state of the vehicle, for example, primarily the roll angle [49], since its dynamics are highly affected by these variables. Therefore, any control algorithm must consider the two-wheeled vehicle roll angle and should be able to estimate this variable using a dynamic filter in case this geometric information is not directly available from the measurement system [39]. More importantly, linear models hold for only small ranges of this variable, in contrast to nonlinear ones. Thus, nonlinear analytic models may offer significant insights in understanding complex dynamics behaviors and develop an efficient model suitable for real-time control outside of the linear regime [50]. On the other hand, in general, the nonlinear dynamic behavior of two-wheeled vehicles is complex. Thus, it is indispensable for model-based control development to carry out a trade-off between model accuracy and its simplicity or usability [47]. For instance, the advanced multibody models can be too large and complex and cannot be used for real-time control purposes except with significant simplifications. Besides, given the direct correlation between model complexity and the vibration modes it can identify [51], clear criteria must be set to determine the level of detail required for a useful dynamical model [32]. Therefore, the derivation of a mathematical model of a two-wheeled vehicle with moderate complexity, which can correctly capture the essential dynamics issues mentioned before, is desirable [30], and these requirements represent the main motivations for the development of the work reported in this paper.

1.3. Literature Review

The development of a good model for a two-wheeled vehicle, as described previously, is not a trivial task indeed. In particular, considering that the predominant linearized equations from the literature are not based on a systematic linearization of full nonlinear differential equations, this task is even more challenging [52,53]. Thus far, systematic linearizations have not achieved analytical expressions for the linearized equation coefficients, until recently, when some authors, such as in [15,20,21], have currently achieved it by developing an automatic computer-aided linearization procedure, finding agreement with the well-known benchmark models [45]. Therefore, due to the traditional approaches found in the literature that were followed during the time, most of the models available are constructed by ad hoc linearization having stability analysis scopes, thereby featuring lots of simplifications, or are based on nonlinear advanced multibody dynamics models for accurate numerical simulations that turn out to be too complicated to handle. This led to a lack of discussion on the relations between the steering and roll angles and the other system key parameters. In particular, their relation with the effective steering and camber angles of the front fork and the contact points between the wheels and the ground, fundamental for the vehicle stability [2], should be further explored. Consequently, a detailed model of a two-wheeled vehicle is complex because this system has many degrees of freedom, and its geometry is intricate [54]. In particular, the mathematical modeling presents considerable difficulties in expressing the position of the front and rear frames [55]. This arises from the front contact geometry, being one of the most complex aspects of the system. In general, establishing the geometric constraint equations is the first challenging problem to study two-wheeled vehicle dynamics [56]. For instance, it is well known that the algebraic equations of these geometric constraints are complex implicit expressions with high nonlinearity [57].

Consequently to the above, a comprehensive, explicit kinematic model of a two-wheeled vehicle with its steering axis set in a general position is not possible, in principle, without resorting to approximations. This is also mentioned by Kane [58], who formulates fundamental kinematic expressions and then resorts to a minimum coordinate approach by substituting explicit linearized expressions. A similar methodology is presented in [2], where the kinematic analysis of a bicycle takes place. The main features of the bicycle position and orientation are expressed as explicit functions of the steering and roll angles. However, this model also uses partial linearization of the variables, being subjected to the same limitations as other models. Both works of Kane and Huang presented a discussion on the formulation of the contact constraint [2,58], particularly on the definition of the position vector of the contact points on the wheels. In contrast, in Frosali and Ricci [59], and in Cossalter [60], advanced nonlinear kinematic models are presented, which proved to be in good agreement with experimental results. However, a discussion on the details of the formulation of these kinematic models and their implication on the vector quantities of interest for the geometric analysis is not explicitly provided. This paper, therefore, tries to address the important issues mentioned before by proposing an advanced kinematic model of the closed-chain geometry of bicycle systems and discuss all the implications of the proposed nonlinear model.

In general, the position of the contact points in a local reference frame attached to the vehicle changes as a function of the system configuration, and this geometric issue is well-established [61]. Therefore, a contact model is required in the kinematic analysis of two-wheeled vehicles since the definition of this model has a significant impact on the numerical results generated by the entire model. However, in the literature, a thorough discussion on this topic is almost systematically eluded [62]. Thus, geometric considerations on the vehicle kinematics are among the main factors complicating the analysis of two-wheeled vehicles. In particular, some rather subtle issues must be addressed, explicitly or implicitly, during the modeling process. Namely, at least two important problems must be pointed out and addressed, that is, (a) the identification of the instantaneous contact points of the wheels given a non zero steering and roll angle, and (b) the formulation of the relation between the steering angle, the roll angle, and rear frame pitch angle [58,63]. Given the problem complexity, most researchers introduced simplifications, like the assumption of proportionality between the steering angle and the roll angle [64], or the assumption that the contact points on the wheels lie in the plane of symmetry of the vehicle [65]. In particular, the latter assumption is also due to the linear approximation, since the front wheel contact point with the road has a constant set of coordinates in the local reference frame of the frontal fork for a linearized model [45]. This paper, on the other hand, tries to relax as much as possible all the simplifications mentioned before to develop a fully nonlinear kinematic model based on a closed-chain multibody approach to vehicle kinematics.

1.4. Scope and Contributions of This Study

The scope and the contributions of this research work can be summarized as follows. This paper deals with the development of a general kinematic model for the large family of two-wheeled vehicles that is fully nonlinear, and, therefore, capable of capturing the complex behavior of this class of articulated mechanical systems, without resorting to any of the simplifying assumptions usually considered in the literature on this important topic. To this end, the strategy employed in this work follows a systematic computational approach, typical of the multibody system dynamics, to analytically formulate an intuitive nonlinear kinematic model of a general two-wheeled vehicle by using a closed-chain kinematic formulation. This is mainly done to minimize in a general scenario the dimension of the set of redundant Lagrangian coordinates necessary for the identification of the three-dimensional geometric configuration of the system under examination. Furthermore, the main objective of the present research work is to analyze and discuss the modeling approaches used for the definition of the wheel-road contact constraints, the assumptions behind their formulation, and the quality of the numerical results produced by their computer implementation. This is done by proposing two kinematic models considering the closed-chain multibody approach devised in this work, by comparing the proposed models with those found in the literature that were independently developed by other research groups, and by considering an additional multibody model constructed by the authors in the MATLAB simulation environment using the SIMSCAPE MULTIBODY software.

The analysis of the wheel-road contact conditions addressed in this investigation allows for providing a timely discussion on the relations between the steering and roll angles, and the other key parameters that are kinematically relevant for the family of two-wheeled systems. In particular, as a peculiar feature of the present research work, two methods for the formulation of the contact constraints are proposed without recurring to the linearization strategy commonly found in the literature. Additionally, a detailed discussion of the displacement of the contact point in the local reference frame of the front fork is provided, thereby allowing for clarifying the intricate relationship of this variable with an arbitrary kinematic configuration of the mechanical system of interest. For this purpose, the kinematic behavior of the vehicle front assembly, which is composed of the frontal fork, the handlebars, and the front wheel, is studied, and its spatial orientation is analyzed by implicitly and explicitly formulating a rotation matrix that is a function of the Lagrangian coordinates of the kinematic models proposed in this work. Thus, the present work lays the foundations for the subsequent formulation of more general dynamic models of two-wheeled vehicles with as much complexity as required by any particular application. For instance, by doing so, the fundamental problems previously described in the development of Advanced Driver-Assistance Systems (ADAS), and other active safety features for two-wheeled vehicles, could be addressed in future research investigations. Therefore, in the present work, the principal emphasis is given to a detailed and entirely reproducible description of the mathematical formulation of a general kinematic model of a two-wheeled vehicle. Also, exhaustive numerical experimentation was carried out for validating the proposed models by means of numerical tests and to pave the way for future developments of dynamic models based on the same closed-chain multibody approach. In particular, a detailed comparison is made between the two contact constraint formulation methods proposed in this work, together with the numerical results generated by the principal analytical approaches taken from the literature. Finally, an additional multibody model numerically implemented by the authors by using the SIMSCAPE MULTIBODY software was also employed for thoroughly benchmarking all the developments of this paper.

1.5. Organization of the Manuscript

This paper is organized as follows. In Section 2, the description of the multibody approach used to derive the kinematic model for two-wheeled vehicles of interest for this paper and the explicit formulation of the Euler angles of the vehicle front assembly are reported. Section 3 presents the formulation of the numerical experiments performed considering the closed-chain multibody model developed in this work, together with its final results and the corresponding discussion on them. Finally, Section 4 includes the summary of the paper, the conclusions reached in this study, and some viable directions for future research.

2. Multibody Approach for the Mechanical Modeling of Two-Wheeled Vehicles as Kinematic Chains

In this section, the fundamental steps of the systematic analytical approach based on the multibody system theory that is employed in the paper for the geometric construction of a general-purpose kinematic model of two-wheeled vehicles are described in detail. To this end, the formulation of a closed-chain geometric structure is discussed and its use is demonstrated for the derivation of the kinematic model of interest for this investigation, whereas its peculiar mechanical features are emphasized.

The remaining part of this section is organized according to the following structure. Section 2.1 describes the main components of the multibody model considered in this investigation, while Section 2.2 reports the relevant geometric aspects of the mechanical system under examination. In Section 2.3, the kinematic analysis of the general multibody model constructed for two-wheeled vehicles is carried out, whereas Section 2.4 demonstrates the practical implementation of the geometric concept based on the closed-chain kinematic structure of this family of mechanical systems. Section 2.5 and Section 2.6, being respectively based on the formulation of a cross-product equation and the introduction of an additional non-generalized coordinate, both necessary for the identification of the contact point collocated on the front wheel, describe the two alternative approaches devised in the paper for the preliminary derivation and the subsequent enforcement of the contact constraint conditions between the front wheel and the road plane. Finally, Section 2.7 deals with the kinematics of the steering point, while Section 2.8 focuses on the definition of a proper set of Euler angles suitable for the identification of the three-dimensional orientation of the vehicle front assembly, which is composed of the frontal fork, the handlebars, and the front wheel.

2.1. System Multibody Model

The two-wheeled vehicle of interest for this research work is modeled as a multibody system composed of the following four rigid bodies:

• The rear wheel.
• The rear frame.
• The front fork/handlebars.
• The front wheel.

Three revolute joints keep the system assembled while it moves freely on the horizontal plane. The permanent contact constraint between each wheel and the plane is also considered. The first revolute joint links the rear wheel and the rear frame. Perpendicular to this, a second revolute joint allows the steering of the front fork/handlebars. Finally, a third revolute joint allows the rotation of the front wheel with respect to the fork/handlebars. The horizontal road plane is contained by the $\mathbf{XY}$ axes of the inertial reference frame. Therefore, the multibody system model of the two-wheeled vehicle considered herein has a seven-dimensional configuration space, as discussed in detail below.

2.2. System Geometry

The geometry of the system under consideration is shown in Figure 1.

The geometric model includes the set of standard geometric parameters used to describe two-wheeled vehicles. These are defined for the system configuration with a zero steering angle and considering an upright position on a horizontal plane. This set of parameters is defined as follows:

• The distance between the contact points of the rear and front wheel, known as wheelbase and denoted with $c_w$.
• The distance between the front contact point and the steering column intersection with the road plane, known as trail and denoted with t.
• The rear and front wheels radii respectively denoted with $a_r$ and $a_f$.
• The tilt angle of the steering column, known as the caster angle and denoted with $\lambda$.

Besides, three additional parameters are used to simplify the kinematic analysis. These are defined as functions of the standard geometrical parameters, as respectively shown in Equations (1)–(3).

$$d=a_{f}sin\left(\lambda \right)-tcos\left(\lambda \right)$$

(1)

$$L=c_w+t-a_{r}tan\left(\lambda \right)$$

(2)

$$H=\frac{a_f-a_r-dsin\left(\lambda \right)}{cos\left(\lambda \right)}$$

(3)

The fork offset denoted with d is the minimum distance between the center of the front wheel and the steering column. In addition, the auxiliary parameters L and H can be geometrically interpreted from Figure 1 since they serve to complete the geometric description of the model. This process allows for obtaining a schematization of the system geometry in a simplified skeleton, as shown in Figure 2.

2.3. Kinematic Analysis

With the geometry defined, the reference frames and the coordinates vector are now considered. The horizontal road plane is defined by the $\mathbf{X}$ and $\mathbf{Y}$ axes of the global reference frame $\mathbf{OXYZ}$. Its $\mathbf{Z}$ axis direction is opposed to the gravity. Three additional local reference frames are used to describe the orientation of the multibody system, as shown in Figure 3.

The reference frame ${\mathbf{O}}^1{\mathbf{X}}^1{\mathbf{Y}}^1{\mathbf{Z}}^1$ is fixed on the contact point between the rear wheel and the road plane. The position vector of its origin is defined as

$${\mathbf{P}}_r={\left[\begin{array}{ccc}{X}_r& Y_r& 0\end{array}\right]}^T$$

(4)

where $X_r$ and $Y_r$ are the coordinates of the contact point on the ground considering the $\mathbf{X}$ and $\mathbf{Y}$ axes of the global reference frame. The orientation of the local reference frame ${\mathbf{O}}^1{\mathbf{X}}^1{\mathbf{Y}}^1{\mathbf{Z}}^1$ is given by the transformation matrix obtained by two successive rotations. That is, the yaw of the rear wheel is given by the angle $\theta$, and the roll is associated with the angle $\chi$. Therefore, its orientation is independent of the rotation of the wheel. The transformation matrix is defined as follows

$${\mathbf{A}}_{1,0}={\mathbf{R}}_z\left(\theta \right){\mathbf{R}}_x\left(\chi \right)$$

(5)

where the generic rotation matrix ${\mathbf{A}}_{i,0}$ represents the transformation matrix between the local reference system labeled with i, that is associated with the generic rigid body i, and the global reference system labeled with 0. Subsequently, one can define the reference frame ${\mathbf{O}}^2{\mathbf{X}}^2{\mathbf{Y}}^2{\mathbf{Z}}^2$, whose origin is fixed to the center of the rear wheel and its orientation is given by the yaw angle $\theta$, the roll angle $\chi$, and the pitch angle $\mu$ of the rear frame, having the following transformation matrix

$${\mathbf{A}}_{2,0}={\mathbf{A}}_{1,0}{\mathbf{R}}_y\left(\mu \right)$$

(6)

Then, the origin of the reference frame ${\mathbf{O}}^3{\mathbf{X}}^3{\mathbf{Y}}^3{\mathbf{Z}}^3$ is fixed at the location shown in Figure 3, and its orientation is given by

$${\mathbf{A}}_{3,0}={\mathbf{A}}_{2,0}{\mathbf{R}}_y^T\left(\lambda \right){\mathbf{R}}_z\left(\psi \right)$$

(7)

where the caster angle $\lambda$ was considered in the rotation matrix ${\mathbf{R}}_y^T\left(\lambda \right)$ in order to have positive values for a backward tilt of the steering column. In summary, the vector of Lagrangian coordinates employed to define the spatial configuration of the multibody system that models the two-wheeled vehicle is

$$\mathbf{q}={\left[\begin{array}{cccccccc}{X}_r& Y_r& \theta & \chi & \mu & \psi & {\varphi }_r& {\varphi }_f\end{array}\right]}^T$$

(8)

where $X_r$ and $Y_r$ respectively represent the $\mathbf{X}$ and $\mathbf{Y}$ coordinates of the rear wheel contact point, $\theta$ is the yaw angle of the rear assembly (rear wheel and rear frame), $\chi$ is the roll angle of the rear assembly, $\mu$ is the pitch angle of the rear frame, $\psi$ is the steering angle of the front fork/handlebars, while ${\varphi }_r$ and ${\varphi }_f$ are the angular displacements due to the rotation of the rear and front wheels, respectively.

2.4. System Kinematic Chain

A kinematic model based on a closed chain is proposed to study the two-wheeled vehicle kinematics, as shown in Figure 4.

The kinematic chain can be directly written by using a vector approach as

$${\mathbf{P}}_r+{\mathbf{r}}_1+{\mathbf{r}}_2+{\mathbf{r}}_3+{\mathbf{r}}_4={\mathbf{P}}_f$$

(9)

In the following, the components of the kinematic chain are described in detail. The points ${\mathbf{P}}_r$ and ${\mathbf{P}}_f$ respectively identify the position vectors of the contact points of the rear and front wheel with the road plane, which have their Cartesian components given by the following expressions

$${\mathbf{P}}_r={\left[\begin{array}{ccc}{X}_r& Y_r& 0\end{array}\right]}^T$$

(10)

$${\mathbf{P}}_f={\left[\begin{array}{ccc}{X}_f& Y_f& 0\end{array}\right]}^T$$

(11)

The vector ${\mathbf{r}}_1$, starting from the contact point ${\mathbf{P}}_r$ to the center of the rear wheel, can be expressed in the global coordinates of the reference frame ${\mathbf{O}}^1{\mathbf{X}}^1{\mathbf{Y}}^1{\mathbf{Z}}^1$ as

$${\mathbf{r}}_1={\mathbf{A}}_{1,0}\left[\begin{array}{c}\begin{array}{c}0\hfill \end{array}\\ \begin{array}{c}0\hfill \end{array}\\ a_r\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}{a}_{r}sin\left(\chi \right)sin\left(\theta \right)\hfill \end{array}\\ \begin{array}{c}-a_{r}sin\left(\chi \right)cos\left(\theta \right)\hfill \end{array}\\ a_{r}cos\left(\chi \right)\end{array}\right]$$

(12)

Subsequently, the vector ${\mathbf{r}}_2$, defined in the reference frame ${\mathbf{O}}^2{\mathbf{X}}^2{\mathbf{Y}}^2{\mathbf{Z}}^2$ with its tail located on the rear wheel center and its head on the point ${\mathbf{O}}_3$, in global coordinates is equal to

$$\begin{array}{c}{\mathbf{r}}_2={\mathbf{A}}_{2,0}\left[\begin{array}{c}\begin{array}{c}L-H{\mathrm{s}}_{\lambda }\hfill \end{array}\\ \begin{array}{c}0\hfill \end{array}\\ H{\mathrm{c}}_{\lambda }\end{array}\right]\hfill \\ \\ =\left[\begin{array}{c}\begin{array}{c}\left(L-H{\mathrm{s}}_{\lambda }\right)\left({\mathrm{c}}_{\mu }{\mathrm{c}}_{\theta }-{\mathrm{s}}_{\chi }{\mathrm{s}}_{\mu }{\mathrm{s}}_{\theta }\right)+H{\mathrm{c}}_{\lambda }\left({\mathrm{c}}_{\theta }{\mathrm{s}}_{\mu }+{\mathrm{c}}_{\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right)\hfill \end{array}\\ \begin{array}{c}\left(L-H{\mathrm{s}}_{\lambda }\right)\left({\mathrm{c}}_{\mu }{\mathrm{s}}_{\theta }+{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\mu }\right)+H{\mathrm{c}}_{\lambda }\left({\mathrm{s}}_{\mu }{\mathrm{s}}_{\theta }-{\mathrm{c}}_{\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right)\hfill \end{array}\\ -{\mathrm{c}}_{\chi }\left(L{\mathrm{s}}_{\mu }-H{\mathrm{c}}_{\lambda -\mu }\right)\end{array}\right]\hfill \end{array}$$

(13)

Then, the vector ${\mathbf{r}}_3$ having its tail fixed at the point ${\mathbf{O}}_3$ and its head located at the center of the front wheel, in global coordinates is written as

$${\mathbf{r}}_3={\mathbf{A}}_{3,0}\left[\begin{array}{c}\begin{array}{c}d\hfill \end{array}\\ \begin{array}{c}0\hfill \end{array}\\ 0\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}d{\mathrm{c}}_{\psi }\left({\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right) - d{\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }\hfill \end{array}\\ \begin{array}{c}d{\mathrm{c}}_{\psi }\left({\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right) + d{\mathrm{c}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }\hfill \end{array}\\ d\left({\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }\right)\end{array}\right]$$

(14)

Finally, two methodologies are studied to write the expression for the position vector ${\mathbf{r}}_4$, which has its tail fixed on the geometric center of the front wheel and its head located at the instantaneous contact point of the front wheel denoted with ${\mathbf{P}}_f$. Each of these methodologies yields a different set of kinematic constraint equations. Therefore, two different kinematic models are proposed in the present work, as described in detail below.

2.5. Cross-Product Model

An effective approach for mathematically formulating the geometric vector denoted with ${\mathbf{r}}_4$ consists of determining the direction of the said vector in the global reference frame, as explained in detail in the following derivations. The unit vector normal to the plane of symmetry of the front wheel can be written in global coordinates as

$${\widehat{\mathbf{n}}}_4={\mathbf{A}}_{3,0}\left[\begin{array}{c}\begin{array}{c}0\hfill \end{array}\\ \begin{array}{c}1\hfill \end{array}\\ 0\end{array}\right]=\left[\begin{array}{c}\begin{array}{c}-{\mathrm{s}}_{\psi }\left({\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right) - {\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\theta }\hfill \end{array}\\ \begin{array}{c}{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{c}}_{\theta }-{\mathrm{s}}_{\psi }\left({\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right)\hfill \end{array}\\ {\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }\end{array}\right]$$

(15)

Then, the unit vector tangent to the front wheel at the contact point can be written as

$${\widehat{\mathbf{t}}}_4=\frac{{\widehat{\mathbf{n}}}_4\times \widehat{\mathbf{K}}}{∥{\widehat{\mathbf{n}}}_4\times \widehat{\mathbf{K}}∥}=\left[\begin{array}{c}\begin{array}{c}\frac{{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{c}}_{\theta } - {\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu } + {\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}\hfill \end{array}\\ \begin{array}{c}\frac{{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }{\mathrm{c}}_{\lambda -\mu } + {\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\theta } + {\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}\hfill \end{array}\\ 0\end{array}\right]$$

(16)

where $\widehat{\mathbf{K}}={\left[\begin{array}{ccc}0& 0& 1\end{array}\right]}^T$ is a constant unit vector opposite to the direction of the gravity acceleration and normal to the road plane. It can be proved that the $\mathbf{Z}$ component of the tangent unit vector ${\widehat{\mathbf{t}}}_4$ is equal to zero. This is a consequence of the fact that the vector ${\widehat{\mathbf{t}}}_4$ must be tangent to the $\mathbf{XY}$ plane. Finally, the vector ${\mathbf{r}}_4$ can be written as

$$\begin{array}{c}{\mathbf{r}}_{4,c}=a_f\left({\widehat{\mathbf{n}}}_4\times {\widehat{\mathbf{t}}}_4\right)\hfill \\ \\ =\left[\begin{array}{c}\begin{array}{c}-\frac{a_f\left({\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi } - {\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }\right)\left({\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }{\mathrm{c}}_{\lambda -\mu } + {\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\theta } + {\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }\right)}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}\hfill \end{array}\\ \begin{array}{c}\frac{a_f\left({\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi } - {\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }\right)\left({\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{c}}_{\theta } - {\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu } + {\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }\right)}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}\hfill \end{array}\\ -\frac{a_f\left(-{\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2 + 2{\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\lambda -\mu }^2 - {\mathrm{c}}_{\chi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi } - {\mathrm{c}}_{\psi }^2 + 1\right)}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}\end{array}\right]\hfill \end{array}$$

(17)

where the subscript c in ${\mathbf{r}}_{4,c}$ indicates that the cross-product approach has been used in the geometric determination of this vector. Since ${\widehat{\mathbf{n}}}_4$ and ${\widehat{\mathbf{t}}}_4$ are mutually perpendicular unit vectors, the magnitude of the vector ${\mathbf{r}}_{4,c}$ is imposed equal to the radius of the front wheel. Each of the vectors computed so far is shown in Figure 5, from which the geometric interpretation of the present formulation approach can be visualized.

With the components of the kinematic chain defined as described above, a closed-loop vector equation can be written as follows

$${\mathbf{P}}_f-{\mathbf{P}}_r={\mathbf{r}}_1+{\mathbf{r}}_2+{\mathbf{r}}_3+{\mathbf{r}}_{4,c}$$

(18)

Then, due to the continuous contact with the road plane, this expression must satisfy the following constraint

$${\widehat{\mathbf{K}}}^T\left({\mathbf{P}}_f-{\mathbf{P}}_r\right)=0$$

(19)

where, considering Equations (12)–(14) and (17), the fundamental equation for modeling the contact condition assumes the following form

$$\begin{array}{c}d\left({\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }\right) + {\mathrm{a}}_r{\mathrm{c}}_{\chi }-{\mathrm{c}}_{\chi }{\mathrm{s}}_{\mu }\left(L-H{\mathrm{s}}_{\lambda }\right) + H{\mathrm{c}}_{\chi }{\mathrm{c}}_{\lambda }{\mathrm{c}}_{\mu }\hfill \\ \\ -\frac{{\mathrm{a}}_f\left(-{\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2 + 2{\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\lambda -\mu }^2 - {\mathrm{c}}_{\chi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi } - {\mathrm{c}}_{\psi }^2 + 1\right)}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}=0\hfill \end{array}$$

(20)

Additionally, the $\mathbf{X}$ and $\mathbf{Y}$ components of Equation (18) are given as follows

$$\begin{array}{c}{X}_f-X_r=\left(L-H{\mathrm{s}}_{\lambda }\right)\left({\mathrm{c}}_{\mu }{\mathrm{c}}_{\theta }-{\mathrm{s}}_{\chi }{\mathrm{s}}_{\mu }{\mathrm{s}}_{\theta }\right) + +H{\mathrm{c}}_{\lambda }\left({\mathrm{c}}_{\theta }{\mathrm{s}}_{\mu }+{\mathrm{c}}_{\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right)\hfill \\ \\ + d\left({\mathrm{c}}_{\psi }\left({\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right) - {\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }\right) + {\mathrm{a}}_r{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\hfill \\ \\ - \frac{{\mathrm{a}}_f\left({\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }\right)\left({\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }{\mathrm{c}}_{\lambda -\mu } + {\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\theta } + {\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }\right)}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}\hfill \end{array}$$

(21)

$$\begin{array}{c}{Y}_f-Y_r=\left(L-H{\mathrm{s}}_{\lambda }\right)\left({\mathrm{c}}_{\mu }{\mathrm{s}}_{\theta }+{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\mu }\right) + H{\mathrm{c}}_{\lambda }\left({\mathrm{s}}_{\mu }{\mathrm{s}}_{\theta }-{\mathrm{c}}_{\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right)\hfill \\ \\ + d\left({\mathrm{c}}_{\psi }\left({\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right) + {\mathrm{c}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }\right) - {\mathrm{a}}_r{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\hfill \\ \\ + \frac{{\mathrm{a}}_f\left({\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi } - {\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }\right)\left({\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{c}}_{\theta } - {\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu } + {\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }\right)}{\sqrt{{\mathrm{s}}_{\lambda -\mu }^2{\mathrm{s}}_{\chi }^2{\mathrm{s}}_{\psi }^2 + 2{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi } + {\mathrm{c}}_{\chi }^2{\mathrm{c}}_{\psi }^2 + {\mathrm{s}}_{\psi }^2{\mathrm{c}}_{\lambda -\mu }^2}}\hfill \end{array}$$

(22)

As will be discussed in detail in the presentation of the numerical results, Equations (20)–(22) allow for fully defining the kinematic configuration of the two-wheeled vehicle system for a set of values of its degrees of freedom. In particular, these expressions are useful to study the behavior of the rear frame pitch angle denoted with $\mu$. Also, the current approach leads to the analytical determination of the relative displacement of the contact points of the wheels as functions of the system degrees of freedom.

2.6. Surface Parametrization Model

An alternative approach to formulate a consistent analytical expression for the vector ${\mathbf{r}}_4$ focuses on noting that, when a set of changes in the degrees of freedom takes place, this fundamental vector experiences an angular displacement in the local reference frame ${\mathbf{O}}^3{\mathbf{X}}^3{\mathbf{Y}}^3{\mathbf{Z}}^3$. This phenomenon can be modeled by introducing an additional geometric variable that serves as a non-generalized coordinate, that is, an angular position coordinate of a massless point that plays the role of an auxiliary geometric parameter. Therefore, the vector of Lagrangian coordinates, when employing this approach, assumes the following different form

$${\mathbf{q}}_s={\left[\begin{array}{ccccccccc}{X}_r\hfill & Y_r& \theta & \chi & \mu & \psi & {\varphi }_r& {\varphi }_f& \beta \end{array}\right]}^T$$

(23)

where $\beta$ denotes the non-generalized coordinate which identifies the angular position of the contact point in the local reference frame, as shown in Figure 6a for the system in the initial configuration and in Figure 6b in the case of an arbitrary configuration.

In both cases considered in Figure 6, the angles of rotation of the rear and front wheel, which are respectively denoted with ${\varphi }_r$ and ${\varphi }_f$, are set equal to zero. Thus, only changes in the coordinates $\chi$, $\psi$, and $\mu$ take place in these scenarios. In fact, as it will be shown later in the paper by means of numerical experiments, the angular displacement $\beta$ is not a function of the angle ${\varphi }_f$. The coordinate $\beta$ allows for analytically expressing the geometric vector ${\mathbf{r}}_4$ as follows

$${\mathbf{r}}_{4,s}={\mathbf{A}}_{3,0}{\left[\begin{array}{ccc}-a_{f}cos\left(\beta \right)& 0& -a_{f}sin\left(\beta \right)\end{array}\right]}^T$$

(24)

where the subscript s in the vector term ${\mathbf{r}}_{4,s}$ indicates that the surface parametrization approach was used. Subsequently, the vector tangent to the wheel at the point of contact is obtained by partially differentiating Equation (24) with respect to the non-generalized angular displacement $\beta$. This vector has the form

$$\begin{array}{c}{\mathbf{t}}_4=\frac{\partial {\mathbf{r}}_4}{\partial \beta }\hfill \\ \\ ={\mathrm{a}}_f\left[\begin{array}{c}-{\mathrm{c}}_{\beta }\left({\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\theta }-{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }\right)\hfill \\ -{\mathrm{s}}_{\beta }\left({\mathrm{c}}_{\psi }\left({\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right) - {\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }\right)\hfill \\ \\ \\ - {\mathrm{c}}_{\beta }\left({\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\theta }+{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }\right)\hfill \\ -{\mathrm{s}}_{\beta }\left({\mathrm{c}}_{\psi }\left({\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right) + {\mathrm{c}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }\right)\hfill \\ \\ \\ {\mathrm{c}}_{\beta }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\lambda -\mu }-{\mathrm{s}}_{\beta }\left({\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }\right)\hfill \end{array}\right]\hfill \end{array}$$

(25)

Then, a constraint equation can be written by considering that this vector must be perpendicular to the unit vector normal to the road plane denoted with $\widehat{\mathbf{K}}$. Therefore, one can write:

$${\widehat{\mathbf{K}}}^T{\mathbf{t}}_4=0$$

(26)

By doing so, the following constraint equation is obtained

$$s_{\beta }\left(s_{\chi }{s}_{\psi }+s_{\lambda -\mu }{c}_{\chi }{c}_{\psi }\right) - c_{\beta }{c}_{\chi }{c}_{\lambda -\mu }=0$$

(27)

It is worth mentioning that, given the non-linearity of Equation (27), to obtain an explicit expression for the angle $\beta$ without resorting to approximations, being this the common path in the literature, it is not a trivial task. However, in the present work, the fully nonlinear expression will be considered. When all the components of the kinematic chain are defined, a closed-loop is proposed, as previously done with the cross-product model. Again, the vector equation associated with the closure of the geometric loop can be written as follows

$${\mathbf{P}}_f-{\mathbf{P}}_r={\mathbf{r}}_1+{\mathbf{r}}_2+{\mathbf{r}}_3+{\mathbf{r}}_{4,s}$$

(28)

Thus, the contact condition of the wheels, as formulated in Equation (19), yields the second kinematic constraint of the surface parametrized model given by

$$\begin{array}{c}d\left({\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }\right) + a_r{\mathrm{c}}_{\chi }+{\mathrm{a}}_f{\mathrm{c}}_{\beta }\left({\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }\right)\hfill \\ \\ - {\mathrm{c}}_{\chi }{\mathrm{s}}_{\mu }\left(L-H{\mathrm{s}}_{\lambda }\right) + H{\mathrm{c}}_{\chi }{\mathrm{c}}_{\lambda }{\mathrm{c}}_{\mu }+{\mathrm{a}}_f{\mathrm{c}}_{\chi }{\mathrm{s}}_{\beta }{\mathrm{c}}_{\lambda -\mu }=0\hfill \end{array}$$

(29)

Finally, the $\mathbf{X}$ and $\mathbf{Y}$ components of Equation (28) are equal to

$$\begin{array}{c}{X}_f-X_r=d\left({\mathrm{c}}_{\psi }\left({\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right) - {\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }\right) + H{\mathrm{c}}_{\lambda }\left({\mathrm{c}}_{\theta }{\mathrm{s}}_{\mu }+{\mathrm{c}}_{\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right)\hfill \\ \\ + {\mathrm{a}}_f{\mathrm{c}}_{\beta }\left({\mathrm{c}}_{\psi }\left({\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }+{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\right) - {\mathrm{c}}_{\chi }{\mathrm{s}}_{\psi }{\mathrm{s}}_{\theta }\right) + {\mathrm{a}}_r{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }\hfill \\ \\ - {\mathrm{a}}_f{\mathrm{s}}_{\beta }\left({\mathrm{s}}_{\lambda -\mu }{\mathrm{c}}_{\theta }-{\mathrm{s}}_{\chi }{\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }\right) + \left(L-H{\mathrm{s}}_{\lambda }\right)\left({\mathrm{c}}_{\mu }{\mathrm{c}}_{\theta }-{\mathrm{s}}_{\chi }{\mathrm{s}}_{\mu }{\mathrm{s}}_{\theta }\right)\hfill \end{array}$$

(30)

$$\begin{array}{c}{Y}_f-Y_r=d\left({\mathrm{c}}_{\psi }\left({\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right) + {\mathrm{c}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }\right) + H{\mathrm{c}}_{\lambda }\left({\mathrm{s}}_{\mu }{\mathrm{s}}_{\theta }-{\mathrm{c}}_{\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right)\hfill \\ \\ + {\mathrm{a}}_f{\mathrm{c}}_{\beta }\left({\mathrm{c}}_{\psi }\left({\mathrm{s}}_{\theta }{\mathrm{c}}_{\lambda -\mu }-{\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\right) + {\mathrm{c}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\psi }\right) - {\mathrm{a}}_r{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }\hfill \\ \\ - {\mathrm{a}}_f{\mathrm{s}}_{\beta }\left({\mathrm{s}}_{\lambda -\mu }{\mathrm{s}}_{\theta }+{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{c}}_{\lambda -\mu }\right) + \left(L-H{\mathrm{s}}_{\lambda }\right)\left({\mathrm{c}}_{\mu }{\mathrm{s}}_{\theta }+{\mathrm{s}}_{\chi }{\mathrm{c}}_{\theta }{\mathrm{s}}_{\mu }\right)\hfill \end{array}$$

(31)

The set of Equations (27), (29)–(31) allows for studying the kinematics of the two-wheeled vehicle system in detail. In particular, when consistent values are assigned to the degrees of freedom of the multibody system, the displacement of the contact point of the front wheel, defined relative to that of the rear wheel, can be readily determined.

2.7. Steering Point Kinematics

As shown in Figure 7, the point in the road plane created by the projection of the steering axis is named the steering point.

This is also a geometric parameter of interest for the mechanical design, as it can provide insight into the kinematics of two-wheeled vehicles in general. Specifically, the magnitude obtained from the expression $∥{\mathbf{S}}_P-{\mathbf{P}}_f∥$, known as the mechanical trail, is also studied for its influence on the system stability [66]. On the other hand, the trajectory of the steering point is studied here, whose position vector is written as a function of the system generalized coordinates as follows

$${\mathbf{S}}_P={\mathbf{P}}_r+{\mathbf{r}}_1+{\mathbf{r}}_2+{\mathbf{r}}_{sp}={\mathbf{P}}_r+{\mathbf{r}}_1+{\mathbf{r}}_2+{\mathbf{A}}_{3,0}\left[\begin{array}{c}0\\ 0\\ -h\end{array}\right]$$

(32)

where the term h in Equation (32) represents the distance from the point ${\mathbf{O}}^3$ to the steering point expressed in local coordinates. This is found by noting that the $\mathbf{Z}$ component of the steering point in global coordinates must be equal to zero. Thus, this consideration yields the following expression

$$h=\frac{{\mathrm{a}}_r-Lsin\left(\mu \right) + Hcos\left(\lambda -\mu \right)}{cos\left(\lambda -\mu \right)}$$

(33)

By back-substituting Equation (33) into Equation (32), the steering point position vector can be written as

$${\mathbf{S}}_P=\left[\begin{array}{c}{X}_r+\frac{cos\left(\theta \right)\left(Lcos\left(\lambda \right) + {\mathrm{a}}_{r}sin\left(\lambda - \mu \right)\right)}{cos\left(\lambda - \mu \right)}\\ Y_r+\frac{sin\left(\theta \right)\left(Lcos\left(\lambda \right) + {\mathrm{a}}_{r}sin\left(\lambda - \mu \right)\right)}{cos\left(\lambda - \mu \right)}\\ 0\end{array}\right]$$

(34)

From this expression, it can be observed that, given a null yaw angle denoted with $\theta$ and when arbitrarily changing the steering angle denoted with $\psi$, the displacement of the steering point is not zero only in the $\mathbf{X}$ component of the system global reference frame.

2.8. Front Assembly Orientation

Due to the method employed to formulate the transformation matrix of the front assembly, the set of Euler angles describing the orientation of such subsystem are not explicitly known. However, these angles can be obtained from the components of the matrix presented in Equation (7) by noting that three successive rotations can fully describe the orientation of the front frame assembly. In particular, this can be done by considering a set of Euler angles based on the sequence Z-X-Y as follows

$${\mathsf{Ã}}_{3,0}=\left[\begin{array}{ccc}{\mathrm{c}}_{\varphi }{\mathrm{c}}_{\delta }{\mathrm{s}}_{\alpha }{\mathrm{s}}_{\varphi }{\mathrm{s}}_{\delta }& -{\mathrm{c}}_{\alpha }{\mathrm{s}}_{\delta }& {\mathrm{c}}_{\delta }{\mathrm{s}}_{\varphi }+{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\alpha }{\mathrm{s}}_{\delta }\\ {\mathrm{c}}_{\varphi }{\mathrm{s}}_{\delta }+{\mathrm{s}}_{\alpha }{\mathrm{c}}_{\delta }{\mathrm{s}}_{\varphi }& {\mathrm{c}}_{\alpha }{\mathrm{c}}_{\delta }& {\mathrm{s}}_{\varphi }{\mathrm{s}}_{\delta }-{\mathrm{c}}_{\varphi }{\mathrm{s}}_{\alpha }{\mathrm{c}}_{\delta }\\ -{\mathrm{c}}_{\alpha }{\mathrm{s}}_{\varphi }& {\mathrm{s}}_{\alpha }& {\mathrm{c}}_{\alpha }{\mathrm{c}}_{\varphi }\end{array}\right]$$

(35)

where $\delta$, $\alpha$, and $\varphi$ represent the absolute yaw, roll, and pitch angles of the front frame, respectively. The explicit expressions for these geometric quantities, as functions of the generalized coordinates, can be obtained by equating the components of Equation (Section 2.8) with those of the matrix reported in Equation (7). For instance, the angle $\alpha$ can be obtained as

$${\mathbf{A}}_{3,0}\left(3,2\right)={\mathsf{Ã}}_{3,0}\left(3,2\right)$$

(36)

which leads to

$$sin\left(\alpha \right) = cos\left(\psi \right)sin\left(\chi \right) - sin\left(\lambda -\mu \right)cos\left(\chi \right)sin\left(\psi \right)$$

(37)

An explicit expression for the camber angle as a function of the generalized coordinates results valuable when considering tire models for two-wheeled vehicles [67,68]. In particular, for computing the tire lateral force, which strongly depends on this parameter, the determination of the camber angle represents a fundamental geometric task which must be fully addressed preliminarily to the dynamic analysis. Similarly, the equation for computing the yaw angle of the front assembly denoted with $\delta$ can be formulated as follows

$${\mathbf{A}}_{3,0}\left(1,2\right)={\mathsf{Ã}}_{3,0}\left(1,2\right)$$

(38)

which leads to

$$\begin{array}{c}sin\left(\delta \right)=\frac{sin\left(\psi \right)\left(cos\left(\theta \right)cos\left(\lambda - \mu \right) + sin\left(\lambda - \mu \right)sin\left(\chi \right)sin\left(\theta \right)\right)}{cos\left(\alpha \right)}\hfill \\ + \frac{cos\left(\chi \right)cos\left(\psi \right)sin\left(\theta \right)}{cos\left(\alpha \right)}\hfill \end{array}$$

(39)

and

$${\mathbf{A}}_{3,0}\left(2,2\right)={\mathsf{Ã}}_{3,0}\left(2,2\right)$$

(40)

which leads to

$$\begin{array}{c}cos\left(\delta \right)=-\frac{sin\left(\psi \right)\left(sin\left(\theta \right)cos\left(\lambda - \mu \right) - sin\left(\lambda - \mu \right)sin\left(\chi \right)cos\left(\theta \right)\right)}{cos\left(\alpha \right)}\hfill \\ + \frac{cos\left(\chi \right)cos\left(\psi \right)cos\left(\theta \right)}{cos\left(\alpha \right)}\hfill \end{array}$$

(41)

Finally, the pitch angle of the front assembly denoted with $\varphi$ can be expressed as a function of the system generalized coordinates as follows

$${\mathbf{A}}_{3,0}\left(3,1\right)={\mathsf{Ã}}_{3,0}\left(3,1\right)$$

(42)

which leads to

$$sin\left(\varphi \right)=-\frac{sin\left(\chi \right)sin\left(\psi \right) + sin\left(\lambda -\mu \right)cos\left(\chi \right)cos\left(\psi \right)}{cos\left(\alpha \right)}$$

(43)

and

$${\mathbf{A}}_{3,0}\left(3,3\right)={\mathsf{Ã}}_{3,0}\left(3,3\right)$$

(44)

which leads to

$$cos\left(\varphi \right)=\frac{cos\left(\chi \right)cos\left(\lambda -\mu \right)}{cos\left(\alpha \right)}$$

(45)

Equations (37), (39), (41), (43) and (45) allow for explicitly defining the orientation of the front assembly of the complete multibody system and can be useful for further modeling tasks in the kinematic and dynamic analysis of two-wheeled vehicles.

3. Numerical Results and Discussion

In this section, a set of numerical experiments is proposed to study the kinematic models developed in this investigation, namely, the two-wheeled vehicle models arising from the use of the cross-product method and the surface parametrization approach. Besides, three additional two-wheeled vehicle models are considered to confront the numerical results found. In particular, two of these kinematic models used for making comparisons are found the literature [59,60], while the third model represents a detailed multibody model created by the authors using the SIMSCAPE MULTIBODY software implemented in the MATLAB computational environment.

The remainder part of this section is organized as follows. Section 3.1 provides a detailed description of the numerical experiments proposed, as well as the numerical values of the physical parameters employed in the two-wheeled vehicle models. Section 3.2 reports the numerical experiment results and the comparative metrics calculated with respect to the model results obtained considering the two-wheeled vehicle model developed using the SIMSCAPE MULTIBODY software. Finally, Section 3.3 provides some final remarks and a detailed discussion on the numerical results found.

3.1. Description of the Numerical Experiments

The purpose and scope of the numerical experiments proposed in this section are to study the system kinematic configuration, employ a kinematically driven approach, and compare the numerical results obtained from this study with those calculated using the kinematic models available in the literature. To this end, in order to implement a kinematically driven approach, the first step consists of analyzing the behavior of the system redundant coordinates when the set of the known values of the degrees of freedom are changed. This is equivalent to analyze the system kinematic constraints, being, in general, these nonlinear algebraic equations among the most complex features in modeling two-wheeled vehicles. In particular, the main goal is to study and understand the relationship between the roll angles of the steering and rear frames with the remaining kinematic variables, primarily resulting from the manipulation of the redundant coordinates in correspondence of the enforcement of the contact between the vehicle wheels and the road. That is, one fundamental objective is to derive the kinematic relationship between the pitch angle $\mu$ and the migration angle $\beta$ associated with the contact point of the front wheel. Similarly, it is of interest for this research to analyze the behavior of the Euler angles of the front assembly when a consistent change in the values of the degrees of freedom is performed. This last kinematic problem does not represent a straightforward task, but it can be readily solved by employing the explicit kinematic expressions found in the present work.

Two-wheeled vehicles are underactuated multibody mechanical systems, which means that the number of external control inputs is less than the number of the system degrees of freedom. Consequently, their complex dynamic behavior cannot be fully described only with a kinematic model, requiring further development of a suitable dynamic model. Besides, if the no-slip condition between the wheels and the ground is considered, although the occupiable regions of the position space remain unaffected, these systems have four nonholonomic constraints and three degrees of freedom in the velocity space [69]. However, the accessible configuration space of the model is still seven-dimensional, as mentioned before. Therefore, to study this family of systems with a kinematically driven approach, seven independent coordinates must be known. For this purpose, it is defined the vector ${\mathbf{q}}_{min}$ with a chosen set of independent coordinates as follows

$${\mathbf{q}}_{min}={\left[\begin{array}{ccccccc}{X}_r\hfill & Y_r& \theta & \chi & \psi & {\varphi }_r& {\varphi }_f\end{array}\right]}^T$$

(46)

where, as already mentioned in the manuscript, $X_r$ and $Y_r$ represent the Cartesian coordinates of the contact point collocated on the rear wheel, $\theta$ and $\chi$ respectively denote the yaw angle and the roll angle of the rear assembly, composed of the rear wheel and the rear frame, $\psi$ represents the steering angle of the handlebars, while the angular displacements denoted with ${\varphi }_r$ and ${\varphi }_f$ are respectively associated with the rotation of the rear and front wheels about their axes. Thus, for a known vector of generalized coordinates ${\mathbf{q}}_{min}$, the kinematic constraints of the model can be numerically solved, leading to appropriate values of the redundant coordinates. In particular, the kinematically driven analysis employed in this work for the model proposed in Section 2.5, namely the cross-product model, consists of the following fundamental steps:

• Substitute the components of the given vector ${\mathbf{q}}_{min}$ into Equation (20).
• Solve numerically this algebraic equation to find the corresponding value of the angle $\mu$ employing the Newton-Raphson method considering an initial estimation ${\mu }_0=0$ (rad) and an appropriate small tolerance denoted with $ϵ$.
• Input the now known vector of Lagrangian coordinates of Equation (8) into Equations (21) and (22) to compute the relative Cartesian coordinates of the contact point of the front wheel.
• Input the vector of Lagrangian coordinates into the explicit expressions given in Equations (37), (39), (41), (43) and (45) for computing the Euler angles of the front assembly.

On the other hand, the surface parametrization model, defined in Section 2.6, includes two kinematic constraints defined in Equations (27) and (29), respectively. Therefore, its augmented vector of Lagrangian coordinates defined in Equation (23) has an additional non-generalized coordinate denoted with $\beta$ to account for the angular migration of the contact point on the front wheel. Thus, the kinematically driven analysis of the parametrized surface model is summarized as follows

• Substitute the components of the given vector of Lagrangian coordinates ${\mathbf{q}}_{min}$ into Equations (27) and (29).
• Solve the set of nonlinear equations resulting from the previous step using the Newton-Raphson method considering the initial estimations ${\mu }_0=0$ (rad) and ${\beta }_0=\lambda$ (rad), as well as a sufficiently small tolerance denoted with $ϵ$.
• Use the now known Lagrangian coordinate vector ${\mathbf{q}}_s$ of the kinematic model to calculate the relative Cartesian coordinates of the front contact point by substituting it into Equations (30) and (31).

In the set of numerical experiments considered in this work, the numerical results of two additional models found in the literature are included to confront the behavior of the proposed analytical approaches. Besides, an additional kinematic model, developed by the authors based on the SIMSCAPE MULTIBODY toolbox belonging to the MATLAB multi-paradigm programming platform, is used as a reference for comparing and validating all the numerical results found in this work. To this end, a performance index based on the Root-Mean-Squared Error (RMSE) is employed as a quantitative metric useful for making consistent comparisons. This index is computed as follows

$$RMSE=\sqrt{\frac{1}{N}\sum _{i=1}^N{\left(y_i-{\widehat{y}}_i\right)}^2}$$

(47)

where N is the number of data, $y_i$ is the reference value, and ${\widehat{y}}_i$ is the experimental data. Figure 8 shows a graphical visualization of the SIMSCAPE MULTIBODY model employed in this work as a reference model.

Figure 8 includes a zoom in the region of contact of the front wheel with the road plane. This shows the approach used to identify the instantaneous point of contact, measure its migration angle $\beta$, and define its global position vector ${\mathbf{P}}_f$. In general, the software SIMSCAPE MULTIBODY allows for performing measurements of the kinematic variables of the rigid bodies defined in the dynamical simulation. However, there is a software limitation of its measurement capacity related to the existence of contact points between rigid bodies. In fact, it is only possible to carry out measurements in the centers of mass or fixed local reference frames of two bodies in contact. However, no information regarding the location of the instant contact point can be obtained directly. Therefore, as a workaround for these software limitations, a very large number of contact proxies were created along the perimeter of the front wheel within the range $\beta =[-85,85]$ (deg), that is, $N_c=1000$, where $N_c$ denotes the number of contact nodes. By doing so, the instantaneous point of contact is identified by the simultaneous measurement of the normal force of the created proxies. In this way, the migration angle measurement is discretized using a fine grid with multiples of $0.003$ (deg). Regarding the orientation angles of the front assembly, on the other hand, SIMSCAPE MULTIBODY allows the direct measurement of numerical values of the transformation matrices associated with each rigid body. Consequently, by implementing a numerical procedure similar to the one described in Section 2.8, these values can be indirectly calculated. Furthermore, the numerical values of the parameters employed for the numerical experiments are reported in

Table 1. Numerical values of the physical parameters employed in the numerical experiments.

Symbol	Meaning	Value (Units)
$c_w$	Wheelbase	$1.02$ (m)
t	Trail	$0.05$ (m)
$\lambda$	Caster angle	30 (deg)
$a_r$	Rear wheel radius	$0.3$ (m)
$a_f$	Front wheel radius	$0.35$ (m)
$ϵ$	Solver tolerance	$1\times {10}^{-8}$ (m) or (rad)

.

The numerical values of the physical parameters reported in Table 1 are consistently used to evaluate the five kinematic models considered in this investigation, namely the SIMSCAPE MULTIBODY kinematic model, the cross-product kinematic model, the surface parameterization kinematic model, the Cossalter kinematic model [60], and the Frosali and Ricci kinematic model [59]. Furthermore, to carry out a consistent comparison, one should note that there is a discrepancy in the conventions for defining the reference frames employed in the kinematic models proposed in [59,60], and these differences should be carefully accounted for in the confront. For the benefit of the interested reader, the key equations of the kinematic models proposed in [59,60] are reported in Appendix A. For this reason, the kinematically driven analysis previously described is carried out for a full rotation of the handlebars of the two-wheeled vehicle models. Thus, the steering angle is varied between $[0,360]$ (deg). Additionally, two values of the roll angle of the rear frame, namely $[0,-15]$ (deg), are considered. Besides, the coordinates $X_r$, $Y_r$, and $\theta$ are set to zero during the numerical experiments to simplify the comparison. This is possible because these three generalized coordinates are of interest only to define the system position with respect to the inertial reference frame, but do not influence the relative configuration of the front assembly or the constraint of permanent contact of the wheels with the road. These assumptions on the measures derived from the numerical experiments make it possible to simplify the adaptations required in the results due to discrepancies in the definition of the reference frames mentioned above.

3.2. Kinematically Driven Study and Comparative Analysis

In this subsection, the numerical results of the kinematically driven study and the corresponding comparative analysis developed in this work are presented. More specifically, the first part of this subsection is devoted to the comparative analysis of the trajectory followed by the contact point belonging to the front wheel, while the second part of this subsection focuses on the study of the set of Euler angles representing the orientation of the front assembly.

The trajectory described by the point of contact of the front wheel on the road plane, given a full rotation of the the handlebars, is shown in Figure 9.

For a zero roll angle, it is observed that the trajectory is symmetric about the $\mathbf{X}$ axis, as shown in Figure 9a. This not being the case when the roll angle is equal to −15 (deg), as evidenced in Figure 9b. In particular, restricting the plotting of the numerical results within the range of the steering angle $\psi =[-60,60]$ (deg) produces the graphical representation depicted in Figure 10.

Within this interval, considering a null roll angle, $21.9\%$ of the total trajectory range is covered on the $\mathbf{X}$ axis, while $11.8\%$ results on the $\mathbf{Y}$ axis. In contrast, having a fixed roll angle equal to $-15$ (deg) increases the covered range to $37.6\%$ and $23.3\%$ in the $\mathbf{X}$ and $\mathbf{Y}$ axis, respectively. In Figure 10a,b, the marker ▷ represents the position of the contact point when the steering angle is zero. Also, in Figure 10a,b, the markers ▽ and indicate the points corresponding to $-60$ (deg) and 60 (deg), respectively.

Given the geometry of the trajectory described by the front wheel contact point, it is of interest to study its similarity with a looped line that is mathematically described by the Limacon curve. Such comparison is presented in Figure 11.

The Limacon curve is described in polar coordinates by the following equation

$$r=a-bcos\left(\psi \right)$$

(48)

where r represents the radial distance of a generic point on the Limacon curve, measured from a reference point set equal to the origin of the global reference system, and $\psi$ denotes the angular displacement from a reference direction associated with the horizontal axis, while the constants a and b define the aspect ratio of the inner loop of the curve. This parametric curve is equivalently described in Cartesian coordinates as follows

$$\mathbf{r}=\left[\begin{array}{c}\xi +cos\left(\psi +\phi \right)\left(a-bcos\left(\psi \right)\right)\\ \eta +sin\left(\psi +\phi \right)\left(a-bcos\left(\psi \right)\right)\end{array}\right]$$

(49)

where $\xi$ and $\eta$ respectively represent the offset in $\mathbf{X}$ and $\mathbf{Y}$ axes, while the parameter $\phi$ accounts for the rotation of the entire curve with respect to the global reference frame. The numerical values of these parameters are found by performing a numerical optimization using the MATLAB built-in function fmincon. To this end, the goal is to minimize the following objective function subjected to a proper set of constraints

$$\left\{\begin{array}{c}min\sum {\left(\widehat{\mathbf{x}}-\mathbf{x}\right)}^2+{\left(\widehat{\mathbf{y}}-\mathbf{y}\right)}^2\hfill \\ \\ \mathrm{subject}\mathrm{to}:\hfill \\ a\ge 0\hfill \\ b\ge 0\hfill \\ a\le b\hfill \end{array}\right.$$

(50)

where $\widehat{\mathbf{x}}$ and $\widehat{\mathbf{y}}$ are the vectors of numerical results of the cross-product model, and $\mathbf{x}$ and $\mathbf{y}$ are the vectors calculated with the use of Equation (49). In addition, the proposed set of algebraic constraints guarantee the existence of the inner loop that characterizes the geometry of the trajectory described by the contact point of the front wheel. For the particular case study considered herein, the optimized numerical values of the parameters are $a=0.1188$, $b=0.1672$, $\phi =1.7370\times {10}^{-10}$, $\xi =1.0676$, and $\eta =7.2167\times {10}^{-11}$. These parameters yield a minimized value of the objective function equal to $1.2650\times {10}^{-4}$ within the interval $\psi =[-60,60]$ (deg).

Figure 12 shows the behavior of the non-generalized parameter $\beta$ when the handlebars is fully rotated.

In the cross-product model, the vector ${\mathbf{r}}_{4,c}$ in Equation (17) is defined a priori in the global coordinate system. Therefore, one does not have an explicit expression for the geometric parameter $\beta$. However, it is found that the third Euler angle of the front assembly, defined in Equations (43) and (45), is closely related to the angular displacement $\beta$. In particular, the numerical values reported in Figure 12 using the cross-product model have been obtained by the following equation

$$tan\left(\beta \right)=-tan\left(\varphi \right)$$

(51)

This expression is obtained from Equations (27), (43) and (45). The results from the numerical approach with the surface parametrized model show a very good agreement. In particular, after considering an offset equal to $\lambda$, the plot in Figure 12 for the surface parametrized model is produced. For the entire rotation of the handlebars, the angle $\beta$ has a range of $69.4912$ (deg) given a null roll angle, and $76.4907$ (deg) when the roll angle is equal to $-15$ (deg). In contrast, when restricting the angular displacement $\psi$ within the interval $[-60,60]$ (deg), the ranges of variation of the geometric parameter $\beta$ are $13.9493$ (deg) for a null roll angle and $28.3808$ (deg) for a roll angle equal $-15$ (deg), respectively. Based on this behavior, the assumption that the front wheel contact point is contained in the symmetry plane of the rear frame, which is equivalent to assuming a null non-generalized coordinate $\beta$, is very frequently considered as a valid hypothesis in the literature. However, for nonlinear analysis, the present results show that this assumption cannot be considered valid in general and represents a non-negligible source of error, as is apparent from the nonlinear kinematic analysis reported and discussed in this work.

Figure 13 presents the behavior of pitch angle of the the rear frame.

It is observed that, for a zero roll angle within the steering interval between 0 (deg) and 45 (deg), the rear frame descends as indicated by the positively increasing trend of $\mu$. Moreover, the range of variation of the angle $\mu$ is equal to $0.1779$ (deg) in this interval. With the steering angle greater than 45 (deg), the rear frame rises, and the angle $\mu$ has negative values. In fact, the plot represented in Figure 13a is symmetric with respect to a vertical line located at the angle $\psi =180$ (deg), where the maximum absolute value of the angle $\mu$ equal to $9.4912$ (deg) is located. In contrast, for a roll angle of $-15$ (deg), as shown in Figure 13b, there is no symmetry. Furthermore, the change in the pitching sense of the rear frame occurs when $\psi =24$ (deg), with a range of variation for the angle $\mu$ in the interval $[0,24]$ (deg) equal to $0.1775$ (deg). In general, the range of variation of the angle $\mu$ is narrow for steering angles of typical magnitudes, i.e., $\psi =[-60,60]$ (deg). This is the reason why linearization approaches for this variable are common in the literature. An expression of this type can be obtained in the framework of the present work by considering Equation (37) to obtain

$$cos\left(\alpha \right)=\sqrt{1-{\left(cos\left(\psi \right)sin\left(\chi \right) - sin\left(\lambda -\mu \right)cos\left(\chi \right)sin\left(\psi \right)\right)}^2}$$

(52)

Then, this expression, which is a function of the Lagrangian coordinates, is substituted into Equations (43) and (45). This yields two explicit expressions for $sin\left(\varphi \right)$ and $cos\left(\varphi \right)$ as follows

$$sin\left(\varphi \right)=-\frac{sin\left(\chi \right)sin\left(\psi \right) + sin\left(\lambda -\mu \right)cos\left(\chi \right)cos\left(\psi \right)}{\sqrt{1-{\left(cos\left(\psi \right)sin\left(\chi \right) - sin\left(\lambda -\mu \right)cos\left(\chi \right)sin\left(\psi \right)\right)}^2}}$$

(53)

$$cos\left(\varphi \right)=\frac{cos\left(\chi \right)cos\left(\lambda -\mu \right)}{\sqrt{1-{\left(cos\left(\psi \right)sin\left(\chi \right) - sin\left(\lambda -\mu \right)cos\left(\chi \right)sin\left(\psi \right)\right)}^2}}$$

(54)

Subsequently, considering the relation reported above in Equation (51), one can eliminate the parameter $\beta$ from Equation (24) by noting that $sin\left(\beta \right)=-sin\left(\varphi \right)$ and $cos\left(\beta \right)=cos\left(\varphi \right)$. In turn, by substituting this analytical result into Equation (29) produces a nonlinear implicit constraint that can be used to solve for the angle $\mu$. Expanding this with first-order Taylor series with respect to the angle $\mu$ and isolating the relevant terms produces

$$\left\{\begin{array}{c}\mu \approx \tilde{\mu }=\frac{d{\sigma }_1 + {\mathrm{c}}_{\chi }\left(a_r + H{\mathrm{c}}_{\lambda }\right) - \frac{a_f}{{\sigma }_2}\left({{\sigma }_1}^2 - {{\mathrm{c}}_{\chi }}^2{{\mathrm{c}}_{\lambda }}^2\right)}{a_f\left({\eta }_2 - {\eta }_1\right) + {\mathrm{c}}_{\chi }\left(\left(L - H{\mathrm{s}}_{\lambda }\right) + d{\mathrm{c}}_{\lambda }{\mathrm{c}}_{\psi }\right)}\hfill \\ \\ {\eta }_1={\mathrm{c}}_{\chi }{\mathrm{c}}_{\lambda }{\sigma }_1\left(\frac{2{\mathrm{c}}_{\psi }}{{\sigma }_2}-\frac{{\mathrm{s}}_{\psi }{\sigma }_1{\sigma }_4}{{\sigma }_3}\right)\hfill \\ \\ {\eta }_2={\mathrm{c}}_{\lambda }{{\mathrm{c}}_{\chi }}^2\left(\frac{2{\mathrm{s}}_{\lambda }}{{\sigma }_2}+\frac{{\mathrm{c}}_{\chi }{{\mathrm{c}}_{\lambda }}^2{\mathrm{s}}_{\psi }{\sigma }_4}{{\sigma }_3}\right)\hfill \\ \\ {\sigma }_1={\mathrm{s}}_{\chi }{\mathrm{s}}_{\psi }+{\mathrm{c}}_{\chi }{\mathrm{c}}_{\psi }{\mathrm{s}}_{\lambda }\hfill \\ \\ {\sigma }_2=\sqrt{1-{{\sigma }_4}^2}\hfill \\ \\ {\sigma }_3={\left(1-{{\sigma }_4}^2\right)}^{3/2}\hfill \\ \\ {\sigma }_4={\mathrm{c}}_{\psi }{\mathrm{s}}_{\chi }-{\mathrm{c}}_{\chi }{\mathrm{s}}_{\lambda }{\mathrm{s}}_{\psi }\hfill \end{array}\right.$$

(55)

where $\tilde{\mu }$ is the approximation obtained for $\mu$. Equation (55) allows for avoiding the numerical iterative solution employed with the use of the cross-product and surface parametrized models. The comparison of the numerical results generated by the previous analytical expression and the SIMSCAPE MULTIBODY model results for the angle $\mu$ is shown in Figure 14.

In general, a good agreement between the proposed approximation and the SIMSCAPE MULTIBODY results is observed in Figure 14. This is confirmed by noting that the RMSE between both results is equal to $0.0257$ for $\chi =0$ (deg), and $0.0420$ for $\chi =-15$ (deg).

The steering point displacement on the $\mathbf{X}$ axis is shown in Figure 15.

A close observation of the numerical results for a zero roll angle and a roll angle equal to $-15$ (deg), represented respectively in Figure 15a,b, reveals a pattern in the effect of that variable. In fact, the presence of a nonzero roll angle eliminates the symmetry concerning a vertical axis located at $\psi =180$ (deg). Moreover, the maximum absolute value of the variable increases. This also holds for the geometric parameters $\beta$ and $\mu$.

Subsequently, Figure 16 represents the numerical results for the system camber angle denoted with $\alpha$, namely the absolute roll angle of the front assembly.

Finally, the remaining two Euler angles of the front assembly are shown in Figure 16 and Figure 17.

In particular, Figure 17 represents the absolute yaw angle of the front assembly denoted with $\delta$. In contrast, the third Euler angle of the front assembly denoted with $\varphi$ is shown in Figure 12 that is based on Equation (51).

As a final remark, in all the numerical results presented in this subsection, a good agreement is found by confronting the proposed analytical approaches based on the cross-product method, the surface parametrization method, the numerical analysis arising from the use of the three-dimensional kinematic model specifically developed in this work using SIMSCAPE MULTIBODY, the independent analytical developments of Cossalter [60], and those by Forsali and Ricci [59] found in the literature. Furthermore, the optimization of the geometric parameters of the Limacon curve served as an additional analytical reference in the representation of the Cartesian coordinates of the contact point located on the front wheel.

3.3. Discussion and Final Remarks

In this subsection, a comprehensive discussion on the numerical results found is provided together with some final remarks. In order to consider a quantitative metric to evaluate and compare the quality of the numerical results, the multibody model developed using the SIMSCAPE MULTIBODY software based on the MATLAB computational environment was considered as the reference model. Assuming the numerical results produced by this multibody model as the reference, the RMSE of the two multibody models proposed in this work, as well as the one of Cossalter [60], and Frosali and Ricci [59], were calculated. The fundamental equations of these two different kinematic models developed by Cossalter, as well as by Frosali and Ricci, are provided in Appendix A. The complete set of numerical values is reported in

Table 2. Calculation of the RMSE for comparing the kinematic models. The SIMSCAPE MULTIBODY simulation results are set as the reference with a fixed roll angle equal to zero.

Model	${\mathit{X}}_{\mathit{f}}$ (m)	${\mathit{Y}}_{\mathit{f}}$ (m)	$\mathit{\beta }$ (deg)	$\mathit{\mu }$ (deg)	$\mathit{\alpha }$ (deg)	$\mathit{\delta }$ (deg)
Cross-prod.	7.0868E−4	$5.8779$E−4	$1.5087$E−1	$1.4063$E−3	$1.3947$E−3	$2.0042$E−1
Sur. param.	$7.0868$E−4	$5.8779$E−4	$1.5087$E−1	$1.4063$E−3	-	-
Cossalter’s	$5.8103$E−3	$8.9220$E−3	$4.0693$	$2.5703$E−2	$5.3286$E−1	$7.1457$E−3
Frosali’s	$7.0293$E−4	$6.0487$E−4	-	$2.0487$E−2	-	-

and

Table 3. Calculation of the RMSE for comparing the kinematic models. The SIMSCAPE MULTIBODY simulation results are set as the reference with a fixed roll angle equal to $-15$ (deg).

Model	${\mathit{X}}_{\mathit{f}}$ (m)	${\mathit{Y}}_{\mathit{f}}$ (m)	$\mathit{\beta }$ (deg)	$\mathit{\mu }$ (deg)	$\mathit{\alpha }$ (deg)	$\mathit{\delta }$ (deg)
Cross-prod.	$7.4645$E−4	$5.6539$E−4	$1.5338$E−1	$1.3769$E−4	$1.3784$E−3	$6.4266$E−4
Sur. param.	$7.4645$E−4	$5.6539$E−4	$1.5338$E−1	$1.3769$E−3	-	-
Cossalter’s	$8.3303$E−3	$1.0530$E−2	$4.3484$	$4.2048$E−2	$6.8937$E−1	$2.3850$E−1
Frosali’s	$7.5011$E−4	$6.1499$E−4	-	$4.2048$E−2	-	-

. In particular, in Table 2, the calculation of the RMSE for comparing the kinematic models is carried out by setting the SIMSCAPE MULTIBODY simulation results as the reference with a fixed roll angle equal to zero, while, in Table 3, the SIMSCAPE MULTIBODY simulation results are obtained with a fixed roll angle equal to $-15$ (deg) for comparing the kinematic models through the computation of the RMSE.

In general, as shown in Table 2 and Table 3, there is good agreement among the five models studied considering the quantitative and qualitative comparison based on their numerical values. In the case of the cross-product and the surface parametrization models, it is pertinent to note that the RMSE values reported in the Table 2 and Table 3 are the same due to the display of the numerical format, which includes only four decimal places. Therefore, for the sake of clarity,

Table 4. Calculation of the RMSE between the cross-product and surface parametrized models for comparing the two models.

$\mathit{\chi }$ (deg)	${\mathit{X}}_{\mathit{f}}$ (m)	${\mathit{Y}}_{\mathit{f}}$ (m)	$\mathit{\beta }$ (deg)	$\mathit{\mu }$ (deg)
0	$3.8706$E−17	$8.4832$E−18	$3.3141$E−15	$1.5307$E−17
15	$3.3009$E−17	$1.3945$E−17	$2.7059$E−15	$7.064$E−17

reports the calculated RMSE between these last two models developed in this work, namely the cross-product and surface parametrized models, for performing a further comparison.

From the analysis of the numerical results presented herein, it is possible to state that the two methods described in this work, namely the cross-product model and the surface parametrization model, are equivalent and were validated through numerical experiments since these were also performed in comparison with other fundamental models found in the literature. However, both the kinematic models developed in this work, being fully nonlinear, do not resort to any simplifying assumptions or linearization procedures, thereby extending their scope of applicability to a wide range of kinematic and dynamic behaviors of the large family of two-wheeled vehicles. Furthermore, the present work presents a detailed description of a systematic and intuitive methodology to formulate kinematic models for two-wheeled vehicles having a general geometric structure. These models are based on the closed-chain kinematic approach, which represents a smart multibody technique capable of avoiding large amounts of redundant coordinates, leading to multibody system simulations that can be easily performed and readily handled considering a small set of redundant coordinates.

Some further comments on the strategy behind the development of the kinematic models are provided below. In general, in the development of the kinematic models presented in this work, the definition of the reference frames turned out to be a fundamental aspect of the model construction since these reference frames were collocated at key geometric points of the two-wheeled system. Besides, the rotation matrices of the reference frames having their origins corresponding to these points were formulated. Then, the analytical expressions for the component vectors of the proposed kinematic chain were derived considering such reference frames. Finally, the wheel contact constraints were formulated, and this step represents another fundamental facet of paramount importance for the kinematic model definition. One particular aspect of the approach employed for the geometric construction followed herein consists of fixing a local reference frame, namely the reference system labeled as ${\mathbf{O}}^1{\mathbf{X}}^1{\mathbf{Y}}^1{\mathbf{Z}}^1$, at the contact point collocated on the rear wheel. Thus, the resulting equations arising from the closed kinematic chain are highly simplified. Furthermore, an advantage of the proposed modeling approach lies in the location of the reference frame ${\mathbf{O}}^3{\mathbf{X}}^3{\mathbf{Y}}^3{\mathbf{Z}}^3$, which allows for formulating the kinematic chain equations in a compact form. This choice proved to be helpful when writing the holonomic constraint of permanent contact between the front wheel and the ground, from which the angle $\mu$ is calculated. In fact, in virtue of such a choice of reference frames, the value of the angle $\mu$ for the system trivial configuration is zero. In general, it is well known that the kinematic loop employed in Equation (19) is highly nonlinear and prevents the reduction of the set of redundant coordinates to a minimal one. Consequently, the numerical solution of the kinematic constraints associated with this geometric feature of the kinematic model was considered in the present study.

At the fundamental intermediate stage of the model definition, two methods were considered to formulate the position vector of the contact point collocated on the front wheel and denoted with ${\mathbf{r}}_4$, giving rise to the two kinematic models presented in this work, namely the cross-product model and the surface parametrization model. These two kinematic models were confronted numerically by performing several experiments and proved to be equivalent, as well as consistent with other kinematic models found in the literature. This being one of the main contributions of this investigation. Additionally, this paper presents an approximation of the geometric angle $\mu$, which is of fundamental importance in vehicle dynamics, by the linearization of the kinematic constraint reported in Equation (29), together with the details of the procedure for its formulation. The resulting expression presents an optimal agreement for large values of the Lagrangian coordinates, as evidenced in Figure 14, as well as by observing the numerical values obtained for the RMSE, namely $0.0257$ for $\chi =0$ (deg), and $0.0420$ for $\chi =-15$ (deg).

Finally, a detailed discussion on the explicit formulation of the position vector of the wheel-ground contact point was presented in this work as well. In particular, the modeling pitfalls concerning the contact point displacement represented in the global and local reference frames were analyzed. The reason for the interest in the study of this fundamental aspect is that it is generally the most complex feature of the kinematic analysis of the tire-road contact mechanics, mainly because a specific material point of the rigid wheel has fixed local coordinates, but the instantaneous point of contact is a function of the Lagrangian coordinates that describe the kinematic model. That is, the position of the instantaneous contact point in a local frame attached to the rear or front frame of the wheels moves when the system configuration changes during the time evolution. However, the simplifying assumption that the local position vector of the contact point will always aim vertically and, consequently, modeling the angular migration of this point as a counter-rotation of the angle ${\varphi }_f$ in the local wheel disk coordinate system, is geometrically incorrect, leading to physically inconsistent numerical results. This simplified approach, which is not used in this work, is common in the development of linearized models, where the hypothesis is acceptable due to the elimination of the higher-order terms in the algebraic equations that describe the kinematic constraints. Therefore, to address this issue, an effective solution that is consistent with the complex geometry of the problem at hand is presented herein, and an analytical expression for the contact point migration angle is systematically derived by noting that it is related to the third angle of the Z-X-Y sequence of Euler angles, specifically defined for this purpose for the front assembly. In particular, for the proposed kinematic model, Equation (51) was found by means of symbolic manipulations and it was subsequently validated numerically through numerical experiments. In this way, it was possible to eliminate the second constraint equation of the surface parametrized model using Equations (53) and (54), thereby proving the independence between the non-generalized parameter $\beta$ and the angle ${\varphi }_f$. To conclude, it is appropriate to mention that the final equation obtained to solve this issue is similar to that proposed by Cosalter [60]. However, the detailed derivation procedure for the solution of this geometric problem is presented herein and its comparison with a nonlinear numerical solution is discussed as well. In fact, as shown Table 2, Table 3 and Table 4, the equation proposed in this work presents a higher agreement with the nonlinear kinematic model developed by the authors in SIMSCAPE MULTIBODY, that was employed as a further external reference. This result, to the best of the authors’ knowledge, was not performed before by other researchers since does not appear in the literature.

4. Summary, Conclusions, and Future Directions of Research

From a general perspective, the authors’ research work is grounded on three main pillars that are multibody dynamics, nonlinear control, and system identification [70,71,72]. In particular, in continuation with the authors’ previous work and for being more specific, this paper treats the fundamental geometric aspects necessary for the development of a general kinematic model, based on a readily viable closed-chain multibody approach, suitable for describing the most relevant kinematic features of a vast class of two-wheeled vehicles.

This paper presents a methodology for the formulation of kinematic models for two-wheeled vehicles, having a general geometry, by using a closed-chain kinematic approach. For this purpose, the most relevant system kinematic features were explored and described through the creation of a proper set of local reference frames that allows for a straightforward derivation of the transformation matrices that define the spatial orientation of the bodies that form the mechanical system, as well as for the geometric construction of its closed-chain geometric shape. Besides, the formulation of global position vectors connecting these frames was performed in such a way so that the geometry of the articulated mechanical system is fully described. In particular, special attention was paid to the methodology employed for deriving the system kinematic constraints modeling the wheel-ground contact, which represents one of the most challenging aspects concerning the mathematical modeling of the closed-chain geometry of two-wheeled vehicles in general. For this purpose, two kinematic models with different methods of formulating the holonomic contact constraint are presented in the paper, that is, an analytical approach based on the cross-product of some fundamental geometric vectors, which appear in the description of the vehicle closed-chain structure, and an analytical technique that leverages the use of a particular non-generalized coordinate, which is specifically introduced to take into account the inherent motion of the contact point on the vehicle front wheel. Additionally, in the last part of the analytical developments presented in this paper, the geometric derivation of the set of Euler angles associated with the vehicle front assembly is performed, and their mathematical form is expressed as an explicit function of the chosen set of Lagrangian coordinates. Furthermore, in the numerical study proposed in the paper, a kinetically driven analysis based on local and global Cartesian and generalized coordinates was proposed to study the kinematic behavior of the contact point collocated on the front wheel. Also, the numerical experiments carried out in the paper, aimed at demonstrating the correctness and effectiveness of the proposed kinematic model, allowed for studying the kinematic behavior of the pitch angle of the vehicle rear frame, as well as the analysis of other kinematic variables of interest. Finally, two paths were followed to establish an appropriate comparison of the proposed kinematic models with the existing ones found in the literature or for confronting other different models specifically devised for this purpose. First, another completely different model was developed by the authors in the MATLAB simulation environment employing the general-purpose software called SIMSCAPE MULTIBODY. Then, two additional kinematic models were taken from the literature for performing meaningful comparisons, namely, those of Cossalter [60], and Frosali and Ricci [59]. These kinematic models, whose mathematical structure is fully available in the literature and is often assumed as a benchmark, were implemented in MATLAB to enable a systematic analysis. By doing so, as shown in details in the numerical results part of the manuscript, an extensive numerical study was carried out in order to validate the proposed kinematic models through numerical experiments, and the Root-Mean-Squared Errors (RMSEs) between the models were calculated to qualitatively and quantitatively confront the numerical results produced by these models.

In summary, the main contributions of this research work include the development of two kinematic models for a general class of two-wheeled systems through the use of a closed-chain kinematic approach that follows the modeling strategy typically employed in multibody systems dynamics. In particular, in the development of this work and in the redaction of the manuscript, the authors gave special attention to reporting in detail all the steps for the formulation process of the characteristic equations that describe the kinematic models illustrated herein. This to present a full and easily reproducible work, allowing also for its use in different model-based developments, where other multibody formulation approaches in redundant coordinates are incompatible or impractical. In this respect, one of the possible viable directions for future research is the development of feedforward and feedback control policies for this family of articulated mechanical systems. In particular, the development of new control systems, where the emerging ADAS applications for two-wheeled vehicles demand nonlinear models with low dimensionality, is one of the most interesting paths for future research [73,74,75]. Additionally, this work provides a discussion confronting two alternative methods suitable for the formulation of the contact constraints. As far as the wheel-road contact is concerned, a campaign of numerical experiments allowed for comparing the numerical results arising from the proposed approach with those obtained from two models taken from the literature and one model specifically developed for this purpose by the authors in MATLAB using the software called SIMSCAPE MULTIBODY. This thorough comparison yielded a numerical verification of the equivalence between the two analytical approaches devised in this investigation, namely the cross-product method and the surface parameterization method. Besides, a detailed discussion on the definition of the local and global position vectors of the contact point on the front wheel was presented in the paper. In particular, the equivalence of the migration angle, associated with the contact point expressed using a set of local coordinates, with the third Euler angle of the Z-X-Y sequence was found. This is a timely discussion, since, during recent years, the problem of the formulation of the holonomic contact constraint for bicycles and motorcycles, despite being non-trivial and giving rise to much of the general difficulties in modeling two-wheeled vehicles, has been systematically eluded in the literature. This is mainly attributable to the literature approach employed during the early stages of the research on two-wheeled vehicles, where ad-hoc linearization approaches were the norm since the full emphasis on the self-stability analysis by eigenvalue criteria was the main topic of interest. This approach was subsequently replaced by multi-degree-of-freedom models having different degrees of complexity, developed by more sophisticated computational tools employing multibody software based on redundant coordinate formulations, and considering contact models that allowed for the separation of the vehicle wheels from the ground with different numerical techniques, such as the penalty method. This transition substantially led to a lack of discussion on the kinematic analysis of the wheel-ground contact suitable for modeling two-wheeled vehicles. Therefore, this issue, which has also recently been highlighted by other authors [2,59,62], was addressed in detail in the present study.

Future works will include further refining the model presented herein, leading to the development of a multibody dynamic model that shares the same closed-chain kinematic approach devised in this investigation. Additionally, once a comprehensive but low-dimensional multibody model is constructed for effectively and efficiently simulating the dynamic behavior of the class of mechanical systems of interest for this study, a systematic formulation of the system equations of motion will be carried out for their subsequent use as a plant in the iterative development and the computer implementation of data-driven nonlinear control strategies, such as those based on the set of Deep Reinforcement Learning (DRL) algorithms [61,70], where the computational cost of a high-dimensional multibody model is restrictive.