The Longitudinal Dynamics of a Road Vehicle: A New Multibody Approach for the Equations of the Rolling Wheel with Constraints Relaxation and Traction Reaction Saturation

Bukoko Ikoki

KDYNAMICS Logiciel, Inc., Sherbrooke, Canada.

DOI: 10.4236/oalib.1109420                

Abstract

The Udwadia-Kalaba formulation is proposed to model the longitudinal dynamics of a road vehicle. To render complex situations such as spinning on a slippery road, an original approach is implemented by the relaxation of constraints in the Udwadia-Kalaba formulation for the rolling of a wheel. In a combined approach of both slip and stiction in the contact section, the constraints equations of pure rolling are associated with stiction. Such constraints are lifted as slip occurs to allow the dynamics of the wheel to take over the normally imposed kinematic contraints. The relaxation of constraints is achieved by the extension of the Udwadia-Kalaba formulation with the semi-least-squares solutions of the constraints equations. This sets biases on the constraints equations based on the description of weight functions that take into account a friction conditionality without branching, which leads to the smooth activation or deactivation of selected constraints equations and associated forces without the need to rewrite the equations of motion.

Keywords

Slip, Stiction, Friction, Constraints, Relaxation, Rolling, Udwadia

1. Introduction

In [1], Udwadia raised the unresolved problem of a proper modeling of the rolling wheel in vehicle dynamics. An attempt [2] is made in this work to address the issue. The complex wheel-and-ground interaction is modeled by considering the rolling constraint. Rolling without slipping is associated with stiction [3]. This constraint needs to be lifted as friction kicks in since a classic variational approach cannot handle dry friction. To a certain extent, the Udwadia-Kalaba formulation for non-ideal constraints [4] takes friction into consideration.

However, by considering the general solution of the least-squares problem associated with constraints equations, it fails to preserve the necessary minimum norm. According to Gauss’ least-constraint principle, the accelerations that are found in this way are not the actual ones. As they remain bound to the constraints equations that they persistently fulfill in a least-square sense, it appears that they do not obey the dynamics of the rolling. Nevertheless, the matrix based U-K formulation appeals in its ability to not only allow the interpretation of the physics behind the equations, conversely, it allows the transposition of an idea into equations.

In [2], we have suggested a way to relax ideal constraints whenever the additional term of non-ideal constraints comes into play in the Udwadia-Kalaba formulation. To allow the relaxation of constraints as slip occurs, the weighted semi-least-squares solutions of the constraints equations are considered instead of the non-weighted constraints equations in the classic formulation. With friction aware weights on the constraints equations, the residuals are altered accordingly. The bias on constraints is achieved to the extent friction which is important by diverting norm minimization efforts to constraints with smaller residuals. By setting the computed generalized coordinates free from selected constraints equations as needed, the true dynamics of the system is allowed to take place.

In their presentation of the joints modeling in flexible multibody systems ( [5], p. 173), Cardona and Géradin address the rolling of the elastic wheel with a slip-stiction approach involving a regularization function proposed by Oden and Martins in ( [6], p. 587). The authors describe such function as problematic with regard to the numerical integration when a perfect zero slipping situation is encountered. However, the retained approach is considered for its simplicity compared to schemes involving constraints activation and deactivation which are regarded as complicated in terms of the time integration procedure. As, according to the authors, they cause violent oscillations of constraints and velocities in the transition phases from the sliding to the stiction cases.

Through the relaxation of constraints, we have presented in [2] a simple method for the simulation of a rolling wheel, which amounts to a smooth activation and deactivation of such constraints, without the inconveniences experienced otherwise. The extended Udwadia-Kalaba equations of motion and the N matrix of constraints relaxation we have introduced, both add a new feature to the multibody dynamics formalism. Such a contribution equally applies to the treatment of a variety of multibody systems that involve intermittent constraints. With our approach, a realistic modeling of ground vehicles can be envisioned. Among other applications, thanks to an accurately modeled target vehicle, an improved automated HIL test procedure for different controllers can be achieved, for both performance and safety goals. While the otherwise involved computations are smoothly performed without upsetting the numerical processes. Thanks to the provided explicit form of the equations of motion of an ODE type which, by an appropriate integration scheme, alleviates the burden of DAE integration and does not require to be rewritten either when constraints vanish. To the best of our knowledge [7], under the Nil novi sub sole provision, no such developments have been presented prior to this work [2].

2. The Equations of Motion of Constrained Systems and Their Relaxation

2.1. The Udwadia-Kalaba Formulation

Udwadia and Kalaba [4] provide the explicit equations of motion of constrained systems by the following formulation:

$\{\begin{array}{l}M\left(q\right)\ddot{q}=Q+Q^c(\text{1a})\hfill \\ A\ddot{q}=b(\; 1b\; )\hfill \end{array}$

M is the n-order positive definite mass matrix, and q is the n dimensional generalized coordinates vector. Matrix A with dimension m × n and m dimensional vector b are obtained from the m constraints functions of m₁ geometric (2a) or m₂ kinematical (2b) kinds.

$\{\begin{array}{l}{h}_g\left(q,t\right)=0(\text{2a})\hfill \\ h_k\left(q,\stackrel{\dot{}}{q},t\right)=0(\; 2b\; )\hfill \end{array}$

Constraint matrix A in (1b) is given by:

$A_{ij}=\{\begin{array}{ll}\partial h_i/\partial q_j,\hfill & i=1,\cdots ,m_1\hfill \\ \partial h_{i-m_1}/\partial {\stackrel{\dot{}}{q}}_j,\hfill & i=m_1+1,\cdots ,m_1+m_2;j=1,\cdots ,n\hfill \end{array}$

Constraints vector b with dimension $m=m_1+m_2$ is obtained from (1b) and (3). Q in (1a) is the matrix of generalized applied forces, conservative or dissipative forces and other complementary inertial forces which do not depend on $\ddot{q}$. Q^c represents the reaction forces that are needed to fulfill the constraints. Constraints (2) are differentiated accordingly to get the second derivative:

$0=\ddot{h}\left(\ddot{q},\stackrel{\dot{}}{q},q,t\right)=A\left(q,t\right)\ddot{q}+\left(\stackrel{\dot{}}{A}\left(\stackrel{\dot{}}{q},q,t\right)+\frac{{\partial }^{2}h\left(q,t\right)}{\partial t\partial q}\right)\stackrel{\dot{}}{q}+\frac{{\partial }^{2}h\left(q,t\right)}{\partial t^2}$ (3)

According to [8], this differentiation amounts to a DAE index reduction that leads to a mild instability of the manifold (3). The Baumgarte technique [9] which would consider unsatisfied (2a) and (2b) as invariant manifolds on the ODE (1a) aims at rendering manifold (3) attracting. By replacing (3) with:

$0=\ddot{h}+2\gamma \stackrel{\dot{}}{h}+{\gamma }^{2}h\text{with}\gamma \ge 0$ (4)

where $\gamma =0$ corresponds to a mildly unstable problem with linearly growing perturbations. Whereas $\gamma >0$ would cause such perturbations to decrease with time. Analytically, the greater the value of $\gamma$, the more attracting the manifold, i.e., asymptotically stable. However, a unanimous grief expressed against Baumgarte is that this property is not verified numerically. But, in the context of this work, the ODE in the UK formulation with the Baumegarte invariants is shown to be part of a stabilization scheme detailed in [10] that effectively addresses this shortcoming for multibody problems with nonholonomic constraints.

2.2. Relaxation of the Constrained Equations

The Gauss least-constraint principle [11] states that the actual motion of a constrained system is the one that minimizes the function:

$G={\left(\ddot{q}-a\right)}^{\text{T}}M\left(\ddot{q}-a\right)$ (5)

subject to some constraints $h\left(q,\stackrel{\dot{}}{q},t\right)=0$. Which after an appropriate number of time differentiations generally result in a linear form identical to (1b). The least-constraint principle clearly amounts to a least-square problem (1b) for a weighted minimum-norm solution ${\Vert \ddot{q}\Vert }_M$ from (5). Building on Gauss’ least-constraint principle, in order to obtain the Udwadia-Kalaba explicit form of the equations of motion [4] subject to the constraints (1b), we set $\ddot{q}=d+a$. Constraints (1b) can then be written as

$A\left(d+a\right)=b\iff Ad=b-Aa\Rightarrow d=X\left(b-Aa\right)+\left(I-XA\right)z$ (6)

Since $\ddot{q}=d+a$ in (6), we have $M\ddot{q}=M\left(d+a\right)=Ma+Md$. Knowing that $a=M^{-1}Q$, we obtain:

$M\ddot{q}=Q+MX\left(b-Aa\right)+M\left(I-XA\right)z$ (7)

Noting by (6) that z is an acceleration, it can conveniently be represented by $M^{-1}C$. Where C is a n dimensional force vector which describes the non-ideal contribution. The constraints are readily obtained in an explicit form in (8a), for the ideal contribution, and (8b) for the non-ideal contribution.

$Q_i^c=MX\left(b-Aa\right)$ (8a)

$Q_{ni}^c=M\left(I-XA\right)M^{-1}C$ (8b)

For the constraints relaxation problem, we suggested in [2] X in the following form:

$X=M^{-\frac{1}{2}}{\left({\left(N^{\frac{1}{2}}\right)}^{+}A{M}^{-\frac{1}{2}}\right)}^{+}{N}^{\frac{1}{2}},N=\text{diag}\left(s_i^2\right)\text{with}i=1,\cdots ,m$ (9)

where N represents the m order positive semidefinite relaxation matrix with weights $0\le s_i\le 1$.

The relaxation of constraints (1b) can be achieved otherwise by considering the more general problem (10) with constraints (10b):

$\{\begin{array}{l}M\left(q\right)\ddot{q}=Q+Q^c(\text{10a})\hfill \\ NA\ddot{q}=Nb(\; 10b\; )\hfill \end{array}$

The constraints forces (8) are reformulated as follows:

$Q_i^c=MY\left(Nb-NAa\right)$ (11a)

$Q_{ni}^c=M\left(I-YNA\right)M^{-1}C$ (11b)

With

$Y=X{N}^{+}=M^{-\frac{1}{2}}{\left({\left(N^{\frac{1}{2}}\right)}^{+}A{M}^{-\frac{1}{2}}\right)}^{+}{\left(N^{\frac{1}{2}}\right)}^{+}$ (12)

where we used the fact that $N^{1/2}{N}^{+}=N^{1/2}{\left(N^{1/2}\right)}^{+}{\left(N^{1/2}\right)}^{+}={\left(N^{1/2}\right)}^{+}$, by virtue of the reverse law applied to generalized inverses for the product of a matrix with its transpose. Matrix Y has the following properties:

$\begin{array}{l}NAYA=NA\\ MYAY=MY\\ {\left(AY\right)}^{\text{T}}N=NAY\\ {\left(YA\right)}^{\text{T}}M=MYA\end{array}$ (13)

For $r\left(A^{\text{T}}NA\right)=r\left(A\right)$, (13) is necessary and sufficient for Y to be A_NM [12]. With {A_NM} being the class of matrices G such that $\hat{x}=Gy$ is a minimum M-seminorm among the N semi-least squares solutions of the possibly inconsistent system $Ax=y$. In other words, $\hat{x}$ minimizes ${\left(y-Ax\right)}^{\text{T}}N\left(y-Ax\right)$ or ${\left(N^{1/2}y-N^{1/2}Ax\right)}^{\text{T}}\left(N^{1/2}y-N^{1/2}Ax\right)$. Also, $\hat{x}=Gy$ is such that ${\Vert y-AGy\Vert }_N\le {\Vert y-Ax\Vert }_N$. This establishes the link with the projection $P_{A\left(N\right)}=AG$ of y onto $\mathcal{R}\left(A\right)$ with respect to ${\Vert y\Vert }_N={\left(y^{\text{T}}Ny\right)}^{\frac{1}{2}}$. As ${\Vert y-P_{A}y\Vert }_N\le {\Vert y-Ax\Vert }_N$. Or ${\Vert N^{1/2}y-N^{1/2}{P}_{A\left(N\right)}y\Vert }_2={\Vert N^{1/2}y-P_{NA\left(I\right)}{N}^{1/2}y\Vert }_2\le {\Vert N^{1/2}y-N^{1/2}Ax\Vert }_2$, $\forall x\in {ℝ}^n$, $\forall y\in {ℝ}^m$, for the system $NAx=Ny$. Knowing that $N^{1/2}{P}_{A\left(N\right)}=P_{NA\left(I\right)}{N}^{1/2}$ according to [12], when N is positive definite, A_NM is the weighted MoorePenrose inverse

$Y=M^{-\frac{1}{2}}{\left(N^{\frac{1}{2}}A{M}^{-\frac{1}{2}}\right)}^{+}{N}^{\frac{1}{2}}$ (14)

which uniquely verifies the four Moore-Penrose conditions:

$AYA=A,YAY=A,{\left(NAY\right)}^{\text{T}}=NAY,{\left(MYA\right)}^{\text{T}}=MYA$ (15)

Provided $r\left(NA\right)=r\left(N\right)$ [13], the constraints (10b) are consistent. Indeed, it can readily be seen that $NAx\in \mathcal{R}\left(NA\right)$. And $Nb\in \mathcal{R}\left(N\right)=\mathcal{R}\left(N\right)A$, because $\mathcal{R}\left(NA\right)\subseteq \mathcal{R}\left(N\right)$ is always true and equality is met when $\dim \left(\mathcal{R}\left(NA\right)\right)=\dim \left(\mathcal{R}\left(N\right)\right)$, meaning $r\left(NA\right)=r\left(N\right)$. Which is assumed here. Consistency is thus established under the provided condition. And because in the context of this work, $A\in {ℝ}_r^{m\times n}$, and $N\in {ℝ}_s^{m\times m}$, where $0<s\le r$, with N diagonal and positive semidefinite, we do have $r\left(NA\right)=r\left(N\right)$ fulfilled. With (10b) consistent, Y provides the M minimum norm solution among all the solutions that actually verify (10b), not in the least-square sense. Y is only required to be an {1, 4M} inverse for which:

$AYA=A,{\left(MYA\right)}^{\text{T}}=MYA$ (16)

Hence, the generalized inverse in (14) would have a simpler expression:

$Y=M^{-\frac{1}{2}}{\left(A{M}^{-\frac{1}{2}}\right)}^{+}$ (17)

Compared with the suggested weighted Moore-Penrose inverses (14) and (17), which are based on assumptions that might well be challenged numerically, it appears (9) has the advantage of encompassing all the situations where N is only required to be positive semidefinite, with possibly inconsistent constraints.

2.3. Stabilization of the Equations of Motion

For the stabilization examination, as a nominal system, we will first consider a system with ideal constraints, with M positive definite and N nonnegative definite. The UK equations can be rephrased in the following state-space form:

$\{\begin{array}{l}\stackrel{\dot{}}{\tilde{q}}=\tilde{v}(\text{18a})\hfill \\ \widetilde{\stackrel{˜}{v}}=\tilde{a}+\tilde{Y}\left(\tilde{b}-\tilde{A}\tilde{a}\right)=\tilde{a}+\tilde{Y}\left(-\frac{\partial \widetilde{\stackrel{˜}{v}}\tilde{v}}{\partial \tilde{q}}\tilde{v}-\tilde{A}\tilde{a}\right)(\; 18b\; )\hfill \end{array}$

where we have set $\tilde{x}=M^{1/2}x$, $\tilde{y}={\left(N^{1/2}\right)}^{+}y$, $\forall x\in {ℝ}^n$, $y\in {ℝ}^m$. And $\tilde{A}={\left(N^{1/2}\right)}^{+}A{M}^{-1/2}$, also $\tilde{Y}=M^{1/2}Y{\left(N^{1/2}\right)}^{+}$. For Y given by (12). By writing

$z=\left[\begin{array}{c}\tilde{q}\\ \tilde{v}\end{array}\right],\hat{h}\left(z\right)=\left[\begin{array}{c}\tilde{g}\left(\tilde{q}\right)\\ \tilde{A}\left(\tilde{q}\right)\tilde{v}\end{array}\right],H={\hat{h}}_z\left[\begin{array}{cc}\tilde{A}& 0\\ \frac{\partial \tilde{A}\tilde{v}}{\partial \tilde{q}}& \tilde{A}\end{array}\right]$ (19)

(18) is stabilized according to [10] by the following scheme

$\stackrel{\dot{}}{z}=\hat{f}-\gamma F\left(z\right)\hat{h}\left(z\right)$ (20)

Which has the same solution as (18) on the manifold $\mathcal{M}\equiv \hat{h}\left(z\right)=0$. Stability [14] is achieved by considering the Lyapunov function $V=\frac{1}{2}{\hat{h}}^{\text{T}}\hat{h}$, and its derivative:

$\stackrel{\dot{}}{V}={\hat{h}}^{\text{T}}\stackrel{\dot{}}{\hat{h}}={\hat{h}}^{\text{T}}H\stackrel{\dot{}}{z}={\hat{h}}^{\text{T}}H\left(\hat{f}-\gamma F\hat{h}\right)$ (21)

By choosing F such that HF is symmetric positive semidefinite with λ₀ as the largest eigenvalue and by determining ${\gamma }_0$ such that $\Vert H\hat{f}\left(z\right)\Vert \le {\gamma }_0\Vert \hat{h}\left(z\right)\Vert ,\forall z$ near $\mathcal{M}\equiv \hat{h}\left(z\right)=0$, exponential stability is found for $\gamma >{\gamma }_0/{\lambda }_0$ as $\stackrel{\dot{}}{V}\le \left({\gamma }_0-\gamma {\lambda }_0\right){\hat{h}}^{\text{T}}\hat{h}<0$. For multibody systems with nonholonomic constraints, we would have (19) as

$\hat{h}\left(z\right)=\left[\begin{array}{c}0\\ \tilde{A}\tilde{v}\end{array}\right],H\hat{f}=\left[\begin{array}{c}0\\ -\tilde{b}+\tilde{A}\tilde{a}+\tilde{A}\tilde{Y}\left(\tilde{b}-\tilde{A}\tilde{a}\right)\end{array}\right]=\left[\begin{array}{l}0\hfill \\ \tilde{c}\hfill \end{array}\right]$ (22)

After development, it can be seen that $\tilde{c}=\left(I-P_{NA\left(I\right)}\right){\left(N^{1/2}\right)}^{+}\left(Aa-b\right)$. Where $P_{NA\left(I\right)}=\left({\left(N^{1/2}\right)}^{+}A{M}^{-1/2}\right){\left({\left(N^{1/2}\right)}^{+}A{M}^{-1/2}\right)}^{+}$. This implies

$\Vert H\hat{f}\Vert \le {\left(N^{1/2}\right)}^{+}\left(Aa-b\right)\le {\lambda }_{\max }\left({\left(N^{1/2}\right)}^{+}\right)\Vert \hat{h}\Vert ={\gamma }_0\Vert \hat{h}\Vert$

By choosing,

$F=\left[\begin{array}{c}0\\ \tilde{Y}\end{array}\right],\text{sothat},HF=\left[\begin{array}{cc}\tilde{A}& 0\\ \frac{\partial \tilde{A}\tilde{v}}{\partial \tilde{q}}& \tilde{A}\end{array}\right]\left[\begin{array}{c}0\\ \tilde{Y}\end{array}\right]=\left[\begin{array}{c}0\\ \tilde{A}\tilde{Y}\end{array}\right]=\left[\begin{array}{c}0\\ P_{NA\left(I\right)}\end{array}\right],$

as $P_{NA\left(I\right)}$ is such that $P_{NA\left(I\right)}=B=B^2=B^{\text{T}}B$ is symmetric and positive semidefinite by definition, exponential stability is established in the sense of Lyapunov for $\gamma >{\gamma }_0/{\lambda }_0$. Baumegarte stabilization is verified by observing that $\stackrel{\dot{}}{z}=\hat{f}-\gamma F\left(z\right)\hat{h}\left(z\right)=\hat{f}-\gamma \tilde{Y}\left(z\right)\hat{h}\left(z\right)$. Which amounts to simply incorporating the stabilizing term $-\gamma \tilde{h}$ into the acceleration compensation factor $\left(\tilde{b}-\tilde{A}\tilde{a}\right)$ in the extended Udwadia-Kalaba fundamental Equation (18b).

The impact of the non-ideal constraints term on stability will be dealt with as a disturbance to the nominal system. For an exponentially stable nominal system,

$\stackrel{\dot{}}{z}=f_n\left(z\right)\text{where}{f}_n=\hat{f}-\gamma F\left(z\right)\hat{h}\left(z\right),$ (23)

with the invariant manifold $\mathcal{M}\equiv \hat{h}\left(z\right)=0$ assimilated to an equilibrium state solution of an ODE, stability of a disturbed system,

$\stackrel{\dot{}}{z}=f_n\left(z\right)+g\left(t,z\right),$ (24)

with $g\left(t,z\right)$ as the disturbing term, is examined here [15] by considering the Lyapunov function $V\left(\hat{h}\right)$ of the nominal system as a candidate in order to verify whether the manifold $\mathcal{M}\equiv \hat{h}\left(z\right)=0$ is an equally stable invariant set for the disturbed system. Stability of the manifod $\mathcal{M}$ is achieved when, for each $ϵ>0$, there is a $\delta \left(ϵ,t_0\right)>0$ such that

$d\left(z\left(0\right),\mathcal{M}\right)<\delta \Rightarrow d\left(z\left(t\right),\mathcal{M}\right)<ϵ,\forall t>t_0$

Because $\stackrel{\dot{}}{V}\le \left({\gamma }_0-\gamma {\lambda }_0\right){\hat{h}}^{\text{T}}\hat{h}=-\lambda V$, the nominal system is exponentially stable with $\hat{h}\left(z\right)=0$ as a stable invariant set. It is readily seen that $V\left(\hat{h}\right)=\frac{1}{2}{\hat{h}}^{\text{T}}\hat{h}$ satisfies ( [15], p. 92)

${\sigma }_1{\Vert \hat{h}\Vert }^2\le V\left(\hat{h}\right)\le {\sigma }_2{\Vert \hat{h}\Vert }^2,\frac{\partial V}{\partial \hat{h}}{f}_n\left(z\right)\le -{\sigma }_3{\Vert \hat{h}\Vert }^2,\frac{\partial V}{\partial \hat{h}}\le {\sigma }_4\Vert \hat{h}\Vert$

Assuming that the perturbation function $g\left(t,z\right)$ has a bound such that

$\Vert g\left(t,z\right)\Vert \le \alpha \Vert \hat{h}\Vert ,\forall t\ge 0\text{with}\alpha \ge 0,$

from $\stackrel{\dot{}}{V}\left(t,\hat{h}\right)=\frac{\partial V\left(\hat{h}\right)}{\partial \hat{h}}{f}_n\left(z\right)+\frac{\partial V\left(\hat{h}\right)}{\partial \hat{h}}g\left(t,z\right)$, we have

$\stackrel{\dot{}}{V}\left(t,\hat{h}\right)\le -{\sigma }_3{\Vert \hat{h}\Vert }^2+\Vert \frac{\partial V}{\partial \hat{h}}\Vert \Vert g\left(t,z\right)\Vert \le -{\sigma }_3{\Vert \hat{h}\Vert }^2+{\sigma }_4\alpha {\Vert \hat{h}\Vert }^2$

Provided $\alpha <\frac{{\sigma }_3}{{\sigma }_4}$, $\stackrel{\dot{}}{V}\left(t,\hat{h}\right)\le -\left({\sigma }_3-{\sigma }_4\alpha \right){\Vert \hat{h}\Vert }^2$ would conclusively establish the manifold $\mathcal{M}\equiv \hat{h}\left(z\right)=0$ as an exponentially stable invariant set according to Khalil approach in [15].

In the context of the problem we are concerned with, we observe that for

$g\left(t,z\right)=\left(I-YNA\right)M^{-1}C$

we have $\Vert \left(I-YNA\right)M^{-1}C\Vert <\Vert M^{-1}C\Vert$, by virtue of the fact that $g\left(t,z\right)$ constitutes the projection of $M^{-1}C$ on the null space of $NA$. Since $M>0$, we also have $\Vert M^{-1}C\Vert \le {\lambda }_{\max }\left(M^{-1}\right)\Vert C\Vert$.

By determining $C_i=\mu N_{ci}\tanh \left(\frac{{\tau }_i}{\mu N_{ci}R}\right)\tanh \left({\hat{h}}_i\right)$ so that $\left|C_i\right|\le \mu N_{ci}\left|\tanh \left({\hat{h}}_i\right)\right|\le {\Gamma }_i\left|{\hat{h}}_i\right|$, with ${\Gamma }_i=\mu N_{ci}$, we finally verify that

$\Vert g\left(t,z\right)\Vert \le \frac{1}{{\lambda }_{\min }\left(M\right)}\Vert \Gamma \Vert \Vert \hat{h}\Vert =\alpha \Vert \hat{h}\Vert$

which qualifies the constraints $\hat{h}=0$ as an exponentially stable invariant set for the disturbed problem with the non-ideal constraints force.

3. The Half Car Model

3.1. Problem Description

A half-car model is represented in Figure 1. The car is to be moved on an occasionally slippery road. A torque is applied progressively on the rear wheel. The front wheel is free to roll. We would like to determine the dynamics of both wheels with respect to the stiction and slip phases.

3.2. The Multibody System

Relative coordinates are used to position the parts. Body 1 is the chassis with three generalized coordinates $q_1,q_2,q_3$. It is characterized by mass $m_1$ and inertia $J_1$. Points $P_0,P_1,P_2$, and $P_3$ are attached to the chassis. Body 2 and 3 are respectively the front arm and the front wheel as shown in Figure 2. The former is positioned with relative generalized coordinate $q_4$, the later with $q_5$. Body 2

Figure 1. Chassis parameters.

Figure 2. Car front arm and wheel parameters.

is characterized by mass $m_2$, inertia $J_2$ and length $L_2$. It bears attachment point $P_4$. Body 3 shows mass $m_3$, inertia $J_3$, stiffness $k_{fw}$, damping $c_{fw}$, and undeformed radius $r_{fw0}$.

Similarly, body 4 and 5 in Figure 3 are respectively the rear arm and wheel. Body 4 with attachment point $P_5$ is positioned by coordinate $q_7$ while $q_6$ positions body 5. The system has 7 degrees of freedom in an open-tree configuration.

In [5], Cardona and Géradin present both the rigid and the flexible body dynamics for the generation of the equations of motion of a multibody system. Figure 4 shows a rigid body. The position $x_P$ of a given point P on a body in the reference frame is expressed here as:

$x_P=x_o+RX$ (25)

X is the vector of the same point P on the body in the local frame. R is the rotation matrix which lines are orthogonal projections of the reference frame base vectors onto the local frame. For rigid bodies in which $\stackrel{\dot{}}{X}=0$, (25) implies

${\stackrel{\dot{}}{x}}_P={\stackrel{\dot{}}{x}}_o+\stackrel{\dot{}}{R}X+R\stackrel{\dot{}}{X}={\stackrel{\dot{}}{x}}_o+\stackrel{\dot{}}{R}{R}^{\text{T}}x={\stackrel{\dot{}}{x}}_o+\tilde{\omega }x={\stackrel{\dot{}}{x}}_o+\tilde{\omega }\left(x_P-x_o\right)$ (26)

Which in the local frame, by premultiplying (26) by $R^{\text{T}}$, corresponds to

$\stackrel{\dot{}}{X}=R^{\text{T}}{\stackrel{\dot{}}{x}}_o+R^{\text{T}}\stackrel{\dot{}}{R}X=R^{\text{T}}{\stackrel{\dot{}}{x}}_o+\tilde{\Omega }X=0\text{where}\tilde{\Omega }=R^{\text{T}}\stackrel{\dot{}}{R}$ (27)

3.3. Constraints Equations

The pure rolling condition is expressed by equating to zero the slipping speed of the wheels. Such a slipping speed is the speed of contact point $x_{Ci}^K$ attached to the road as measured in non-rotating frame of the wheel frame K.

${\stackrel{\dot{}}{x}}_{Ci}^I={\stackrel{\dot{}}{x}}_{Oi}^I+\stackrel{\dot{}}{R}{R}^{\text{T}}{x}_{Ci}^K={\stackrel{\dot{}}{x}}_{Oi}+{\tilde{\omega }}_i\left(x_{Ci}^I-x_{Oi}^I\right)$ (28)

For the rear wheel, ${\stackrel{\dot{}}{x}}_{Oi}={\stackrel{\dot{}}{x}}_{P_5}={\stackrel{\dot{}}{x}}_{P_0}+{\stackrel{\dot{}}{R}}_{ra}{X}_{P_5}$. With $x_{P5}=x_{P_0}+R_{ra}{X}_{P_5}$. For

Figure 3. Car rear arm and wheel parameters.

Figure 4. Body with point P in local frame vs. reference frame.

the front wheel, ${\stackrel{\dot{}}{x}}_{Oi}={\stackrel{\dot{}}{x}}_{P_4}={\stackrel{\dot{}}{x}}_{P_1}+{\stackrel{\dot{}}{R}}_{fa}{X}_{P_4}$, with ${\stackrel{\dot{}}{x}}_{P_1}={\stackrel{\dot{}}{x}}_{P_0}+{\stackrel{\dot{}}{R}}_1{X}_{P_1}$. And $x_{P_4}=x_{P_1}+R_1{X}_{P_4}$. The constraints equations for both wheels are (28):

$\{\begin{array}{l}\stackrel{\dot{}}{h}\left(1\right)\equiv {\stackrel{\dot{}}{x}}_{P4}\left(1\right)-\left({\stackrel{\dot{}}{q}}_3+{\stackrel{\dot{}}{q}}_4+{\stackrel{\dot{}}{q}}_5\right)x_{P4}\left(3\right)=0\text{(29a)}\\ \stackrel{\dot{}}{h}\left(2\right)\equiv {\stackrel{\dot{}}{x}}_{P4}\left(3\right)=0\text{(29b)}\\ \stackrel{\dot{}}{h}\left(3\right)\equiv {\stackrel{\dot{}}{x}}_{P5}\left(1\right)-\left({\stackrel{\dot{}}{q}}_3+{\stackrel{\dot{}}{q}}_6+{\stackrel{\dot{}}{q}}_7\right)x_{P5}\left(3\right)=0\text{(29c)}\\ \stackrel{\dot{}}{h}\left(4\right)\equiv {\stackrel{\dot{}}{x}}_{P5}\left(3\right)=0\text{(29d)}\end{array}$

The velocity-based constraints matrix A is obtained by deriving (2b) once with respect to $\stackrel{\dot{}}{q}$.

$A_{ij}=\frac{\partial {\stackrel{\dot{}}{h}}_i}{\partial {\stackrel{\dot{}}{q}}_j}$ (30)

Assuming scleronomic constraints, vector b in the motion Equation (1) is then obtained from (3) as follows

$b_i=-{\displaystyle {\sum }_{j=1}^n{\displaystyle {\sum }_{k=1}^n\frac{\partial A_{ij}}{\partial q_k}{\stackrel{\dot{}}{q}}_j{\stackrel{\dot{}}{q}}_k}}$ (31)

The constraints-weights matrix N is given by:

$N=\left[\begin{array}{cccc}{s}_f^2& 0& 0& 0\\ 0& s_r^2& 0& 0\\ 0& 0& s_f^2& 0\\ 0& 0& 0& s_r^2\end{array}\right],s_i=1-{\tanh }^2\left(\frac{1}{3}{k}_s\frac{T_{mi}+T_{bi}}{\mu N_{ci}{r}_i}\right),i=\left\{f,r\right\}$ (32)

The constraints equations in (29) are selectively activated, deactivated or simply discriminated against by the stiction coefficients $s_i$. Where $i=\left\{f,r\right\}$ stand respectively for the front and rear wheel. $N_{ci}$ is the contact reaction force of the ground for each wheel. And $r_r=x_{P5}\left(3\right)$ and $r_f=x_{P4}\left(3\right)$, the heights of the wheels center. The matrix C of generalized forces for nonideal constraints is

$C={\left[\begin{array}{ccccccc}{F}_{ci}& 0& 0& 0& 0& 0& 0\end{array}\right]}^{\text{T}}$

$F_{ci}$ represents the sliding friction reaction forces on the front and rear wheels.

$F_{ci}=-\mu N_{ci}\tanh \left(\frac{1}{3}{k}_F\frac{T_{mi}+T_{bi}}{\mu N_{ci}{r}_i}\right)\tanh \left({\omega }_i{r}_i-{\stackrel{\dot{}}{q}}_1\right),i=\left\{f,r\right\}$ (33)

4. Simulations Results

A torque $T_{mot}$ is increasingly applied on the rear wheel. The front wheel is free to roll as the vehicle enters a slippery phase. Subsequently, a braking torque is applied on both wheels. The rotational speeds of both wheels are observed. The applied torques are represented in Figure 5(a) and are given by:

$T_{\text{mot}}=T_{\max }\tanh \left(k_{M}t\right)$ (34)

for the driving torque. And the braking torque is

$T_b=-\frac{T_{bm}}{4}\left(1+\frac{1-{\text{e}}^{-t_{br}\left(t-t_b\right)}}{1+{\text{e}}^{-t_{br}\left(t-t_b\right)}}\right)\left(1+\frac{1-{\text{e}}^{-t_{br}\left(t-t_b-b_d\right)}}{1+{\text{e}}^{-t_{br}\left(t-t_b-b_d\right)}}\right)\tanh \left(\omega \right)$ (35)

The friction coefficient profile, as shown in Figure 5(b) , is

$\mu \left(t\right)={\mu }_l+\left({\mu }_u-{\mu }_l\right)\cdot \tanh {\left(\frac{1}{s_d}\left(t-t_s\right)\right)}^4$ (36)

Air resistance is taken into account with the bearing resistance according to:

$F_w=-\frac{1}{2}\rho A_r{c}_x{\stackrel{\dot{}}{q}}_1^2,T_{be}=-k_{be}{\omega }_i^2,i=\left\{f,r\right\}$ (37)

For a maximum applied torque of $T_m=150\text{N}\cdot \text{m}$, the Kinematic results of both the rear and the front wheels are shown in Figure 6. For the rear wheel, a clear increase in speed is observed as the wheel experiences a much lower friction

(a) (b)

Figure 5. Simulation input. (a) Applied torques; (b) Friction profile.

(a) (b)

Figure 6. Kinematic results for $T_{\max }=150\text{N}\cdot \text{m}$. (a) Rotational speeds; (b) Car speed.

coefficient phase. The traction is displayed in Figure 8(a). The traction reaction $Q_{ci}$ is saturated by the friction limit $\mu N$ according to (33). For a higher maximum applied torque $T_{\max }=550\text{N}\cdot \text{m}$, the kinematic results are shown in Figure 7. A discrepancy is clearly observed in Figure 7(a) at lower rotational speeds between the rear and the front wheel, though the friction coefficient is relatively high. A greater spinning is observed at a lower friction coefficient value. A distinctly saturated reaction force $Q_{ci}$ is displayed in Figure 8(b).

The numerical values for the simulation are presented in Table 1.

A much higher torque might be applied to the rear wheel with $T_{\max }=700\text{N}\cdot \text{m}$. As shown in Figure 9, the traction reaction force $Q_c$ is saturated to the bottom. The sliding friction reaction is saturated by the friction limit $\mu N$.

Table 1. Numerical values.

(a)(b)

Figure 7. Kinematic results for $T_{\max }=550\text{N}\cdot \text{m}$. (a) Rotational speeds; (b) Car speed.

(a)(b)

Figure 8. Traction $Q_i^c$ and slip $Q_{ni}^c$ reactions on the rear wheel. (a) $T_{\max }=150\text{N}\cdot \text{m}$ ; (b) $T_{\max }=550\text{N}\cdot \text{m}$.

The front wheel reaction forces of both traction and sliding friction are obtained from the related torques constraint associated with the rotation of the wheel by division by the wheel radius.

(a)(b)

Figure 9. Simulation results for $T_{\max }=700\text{N}\cdot \text{m}$. (a) Vehicle speed; (b) Reaction forces with limit.

Figure 10 shows the existence of a lower tractive reaction on the front that is saturated by the friction limit when an applied torque on the rear wheel has the maximum value of $T_{\max }=150\text{N}\cdot \text{m}$. Conversely, the friction reaction force is nonexistent since the wheel is rolling freely and is in a stiction mode.

Figure 11 shows the traction reaction force on the front which is initially saturated by the friction limit and subsequently goes down as the vehicle is being slowed down by a more accentuated slip created by a higher torque on the rear wheel with $T_{\max }=550\text{N}\cdot \text{m}$.

Figure 10. Front wheel reaction forces and limit for $T_{\max }=150\text{N}\cdot \text{m}$.

Figure 11. Front wheel reaction forces and limit for $T_{\max }=550\text{N}\cdot \text{m}$.

Figure 12 shows the traction and sliding friction reaction forces on the front wheel for $T_{\max }=700\text{N}\cdot \text{m}$ on the rear wheel. Interestingly, the traction reaction force on the front is not bounded to zero but displays a negative value.

Contrary to the rear wheel which is under an applied torque of the highest maximum value $T_{\max }=700\text{N}\cdot \text{m}$ and is mostly in a slip mode, the front wheel remains in stiction mode because there is no applied torque on it. Thus, as the rear traction reaction is zeroed, the wind force $F_w$ takes over the dynamics of the vehicle because of inertia, and slows it down. This is naturally reflected on the front wheel by a braking effort since it is in a stiction mode. Just as any pushing-back force on the vehicle would generate a negative traction force on a free wheel in a stiction mode. Which would cause the braking and ultimately the reversal of its movement. With a null sliding friction reaction force.

Effect of Baumgarte Parameter Variation on Simulation

Figure 13 shows the effect of increased Baumegarte parameter $\gamma$ on the rotational speeds of both wheels. The reference case is the one considered in Figure 7(a) which uses $\gamma =20$ for a $T_{\max }=550\text{N}\cdot \text{m}$. Figure 13(b) shows the results for $\gamma =100$.

Further increment of Baumegarte parameter to $\gamma =200$ in Figure 14(a) and $\gamma =2000$ as shown in Figure 14(b) displays the gradually attenuated relaxation of constraints and the accentuated compliance outside the relaxation zone. Which might not be desired as it shadows the effect of higher torques in the

Figure 12. Front wheel reaction forces and limit for $T_{\max }=700\text{N}\cdot \text{m}$.

(a)(b)

Figure 13. Vehicle wheel rotational speeds. (a) $\gamma =20$ ; (b) $\gamma =100$.

relaxation of constraints. However, the attractiveness displayed is indicative of the robustness of the integration scheme with a simple choice of the stabilization parameter that could well be calibrated on experimental results.

(a)(b)

Figure 14. Vehicle wheel rotational speeds. (a) $\gamma =200$ ; (b) $\gamma =2000$.

5. Conclusions

In this work, we presented a simple method for the simulation of a rolling wheel for the longitudinal dynamics of a road vehicle.

A new approach for the realistic simulation of the rolling of the elastic wheels is introduced. Such a realistic simulation includes situations such as the spinning of the wheel, which may originate from the encounter with a lower friction zone or the application of a relatively higher torque with respect to the friction limit. These situations are associated with a breakaway from the constraints equations for pure rolling. To account for such a breakaway, an extended formulation of the Udwadia-Kalaba equations of motion for constrained systems was proposed. The method amounts to a smooth activation and deactivation of constraints that are associated with pure rolling according to the encountered situation of stiction or slip. This is done without the inconveniences experienced otherwise in the referenced literature.

The relaxation of a pure rolling constraint naturally leads to the saturation of the associated constraint force which corresponds to the traction force saturation. Thanks to the extended U-K formulation, the computation of the explicit expression of such a traction force is obtained as a result of the dynamics of the rolling wheel, as opposed to the approach that sets it as an input to the model via some approximate formula.

This in our view opens the path to a whole different philosophy in the design of controllers for wheel slip that a future work will explore, where slip in itself is not the direct object for control rather the consequence of an appropriately assessed traction reaction effort with regard to the road condition.

Funding

This research received no external funding.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding author.

Appendix

The simulation code can be obtained on author’s website by scanning the following QR code:

Conflicts of Interest

The author declares no conflicts of interest.

References

[1]	Caverley, R. (2001) Constrained or Unconstrained, That Is the Equation. USC News, September.
[2]	Ikoki, B.K. (2021) The Rolling Wheel Equations without Magic: A Combined Slipstiction Approach by the Udwadia-Kalaba Formulation with Constraints Relaxation for a Smooth Simulation. Journal of Transportation Technologies, 11, 378-389. https://doi.org/10.4236/jtts.2021.113024
[3]	Neimark, J.I. and Fufaev, N.A. (1972) Dynamics of Nonholonomic Systems, Vol. 33 of Translations of Mathematical Monographs. American Mathematical Society, Providence.
[4]	Udwadia, F.E. and Kalaba, R. (1992) A New Perspective on Constrained Motion. Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences, 439, 407-410.
[5]	Cardona, A. and Géradin, M. (2001) Flexible Multibody Dynamics. Wiley, Hoboken.
[6]	Oden, J.T. and Martins, J. (1985) Models and Computational Methods for Dynamic Friction Phenomena. Computer Methods in Applied Mechanics and Engineering, 52, 527-634. https://doi.org/10.1016/0045-7825(85)90009-X
[7]	Soltakhanov, S.K., Yushkov, M.P. and Zeghzda, S. (2009) Mechanics of Nonholonomic Systems. Springer, Berlin. https://doi.org/10.1007/978-3-540-85847-8
[8]	Ascher, U., Honshen, C., Petzold, L. and Sebastian, R. (1995) Stabilization of Constrained Mechanical Systems with DAEs and Invariant Manifolds. Mechanics of Structures and Machines, 23, 135-157. https://doi.org/10.1080/08905459508905232
[9]	Baumgarte, J. (1972) Stabilization of Constraints and Integrals of Motion in Dynamical Systems. Computer Methods in Applied Mechanics and Engineering, 1, 1-16. https://doi.org/10.1016/0045-7825(72)90018-7
[10]	Ascher, U. and Chin, H. (1994) Stabilization of DAEs and Invariant Manifolds. Numerische Mathematik, 67, 131-149. https://doi.org/10.1007/s002110050020
[11]	Kalaba, R., Natsuyama, H.H. and Udwadia, F.E. (2004) An Extension of Gauss’ Principle of Least Constraint. International Journal of General Systems, 33, 63-69. https://doi.org/10.1080/0308107031000139996
[12]	Mitra, S.K. and Rao, C.R. (1974) Projections under Seminorms and Generalized Moore Penrose Inverses. Linear Algebra and Its Applications, 9, 155-167. https://doi.org/10.1016/0024-3795(74)90034-2
[13]	Stanimirovic, P., Pappas, D., Katsikis, V. and Cvetkovic, M. (2015) Outer Inverse Restricted by a Linear System. Linear and Multilinear Algebra, 63, 2461-2493. https://doi.org/10.1080/03081087.2015.1019200
[14]	Lyapunov, A.M. (1907) Problème général de la stabilité du mouvement. Annales de la Faculté des Sciences de Toulouse, 9, 203-474. https://doi.org/10.5802/afst.246
[15]	Khalil, K.H. (2015) Nonlinear Control. Pearson Education Limited, London.