自行车动力学建模及稳定性研究分析

关键词： 匀速 多刚体动力学 稳定性分析 自行车 非完整约束

自行车具有结构简单，使用方便等特点，同时它的动力学行为与一些复杂的力学和数学问题密切相关.静态来看，自行车很难通过前后轮与路面的接触实现稳定直立.因此，其物理本质与“倒立摆”结构类似，都属于一类静态不稳定系统[2].但当自行车在适当的运动速度条件下，却能够实现动态自稳定的直立“自行”.这种自行稳定性现象在历史上引起了众多学者的研究兴趣，并出现了多种多样的理论解释.Rankine[3]将自行车简化为只包含车轮偏转角的单自由度动力学系统，粗略地认为自行车的平衡源于骑行者对车把转角的控制；Bourlet[4]建立了考虑车轮陀螺效应的自行车简单动力学模型，并认为自行车的自稳定性与车轮的陀螺效应有关；文献[5,6,7,8]也建立了各自的自行车等效动力学模型，但这些模型在反映自行车的自稳定特征时均存在一定程度上的不足.1899年法国数学家Carvallo[9]和英国剑桥大学Whipple[10]应用刚体动力学理论，并引入一些线性化的假设，分别独立地建立了自行车做匀速直线运动时的线性化摄动动力学模型.特别是，两位研究者首次将自行车简化为由4个刚体组成的简单构型，并给定了标准的几何和物理参数.该简化的自行车构型被后来的研究者广泛采用[11,12,13,14,15]，称之为“Carvallo-Whipple自行车构型”或简称为“Whipple构型”.

随着对自行车研究的不断深入，以及动力学理论和计算机技术的发展，自行车的研究范畴逐步扩大.当前，关于自行车的研究内容不仅涉及经典力学[16,17,18,19]的多个分支，而且与应用数学[20,21]，控制理论[22,23,24,25]，机器人[26]等多个学科领域密切相关.

论文的结构如下：第2节简要描述了Whipple自行车构型，并对自行车动力学研究中涉及到的完整约束和非完整约束进行了概述；第3节系统总结了自行车动力学建模中的各类方法；第4节对自行车3种动态平衡状态的稳定性问题的研究进行了总结；最后，展望了自行车动力学的研究前景.

1、自行车构型描述及约束方程

图1为通常采用的Carvallo-Whipple自行车构型系统.该系统主要由前后框架H和B(骑行人作为刚体固结于后框架)、前车轮F和后车轮R四个刚体组成.忽略铰接处的间隙、摩擦、阻尼等因素，并在建模过程中考虑以下几个方面的假设：

(1)前、后轮边缘为刀刃型结构，与地面的接触关系为点接触；

(2)前、后轮是圆周对称结构，其质心与各自圆心重合；

(3)假设前轮和后轮与地面之间的最大静摩擦系数足够大，以保证前后轮运动时处于纯滚动状态，即与地面无任何滑移；

(4)前框架和后框架分别关于前轮和后轮所在的平面对称，一是几何对称，二是质量分布对称.

为了描述自行车的运动状态，定义空间惯性坐标系FI={O-ijk}，并将坐标原点设置在后车轮与路面的接触点P处，i轴的正向指向前车轮，k轴正向背离地面，j轴正向指向骑行者的左侧.分别在4个刚体质心处建立固连坐标系

图1Carvallo-Whipple自行车构型图

自行车的4个刚体通过3个理想铰链连接，当不考虑路面对车轮的约束时，描述自行车自由运动需要9个(4×6-3×5=9)广义坐标.Xiong等[27,28]将9个广义坐标定义如下：选择前叉转轴与前框架坐标轴的交点D为参考点，其在惯性坐标系FI中的3个坐标分量(x,y,z)作为系统的3个广义坐标.遵循3-1-2的旋转顺序，选择后框架固连坐标系Fb相对于惯性坐标系FI的3个欧拉角(ψ,θ,?)为系统的另外3个广义坐标，其中ψ为绕着z轴的偏航转角，θ为绕着x轴的倾斜角，?为绕着y轴的俯仰角.另外3个广义坐标分别为车把转角δ(前框架相对于后框架的转角)，后、前轮相对于后、前框架的转角φr和φf.这9个广义坐标可用向量表示为：

另外，描述Whipple自行车构型需9个几何参数和16个质量参数.这9个几何参数包括：自行车的轴距w，前轮尾迹c(即前叉转轴延长线与地面的交点和地面接触点Q之间的距离，当该交点在Q点右侧时，c为正值，否则为负值)，前叉转轴相对于竖直线的倾斜角λ，前后轮的半径(rf,rr)，前框架质心位置(xh,yh)和后框架质心位置(xb,yb).16个质量惯性参数包括：4个刚体的质量(mb,mh,mr,mf)以及它们相对于各自质心的惯量张量(Ib,Ih,Ir,If)(共有4×2+2×2=12个非零分量).Whipple自行车构型的经典参数可参考文献[17].根据车轮纯滚动的假设条件，容易理解自行车是一个同时受完整和非完整约束的多刚体系统.然而，与Chaplygin雪橇[29]和独轮机器人[30]等传统的非完整系统相比，自行车中的约束方程不易用显式方程来描述.由于通过铰链约束连接的前后车轮的轮心和前后框架之间形成了复杂的链式传递关系，接触点处的几何约束往往是隐式和非线性的.为建立Whipple自行车构型下的线性化动力学模型，Meijaard[17]对非线性约束方程进行了线性化近似处理；Psiaki[13]利用解析几何的方法首次推导了完整约束方程的解析形式，其为四次代数方程，且具有高度非线性；Basu-Mandal等[18]基于寻找两个车轮的最低点为参考条件推导了约束方程；类似地，Peterson和Hubbard[31,32]考虑轮胎的形状，采用直观的解析几何方法建立了完整和非完整约束方程；Wang等[33]对两个车轮圆周进行参数化，利用极值条件(车轮与地面的接触点为轮缘的最低点)推导出车轮的完整约束方程.

理论上，接触点处的几何约束方程与两个接触体的几何外形密切相关.因此，在建立由点接触引起的几何约束方程时，需要对两个接触体的几何外形进行参数化处理，进而依赖接触条件，确定接触点处的外形参数与系统广义坐标之间的函数关系，在此基础上，建立由系统广义坐标表示的几何约束方程[30,34].基于Zhao和Liu[34]关于接触约束的理论描述，Xiong等[27,28]分别建立了自行车在水平面和旋转曲面上的几何约束方程和速度约束方程.当自行车在水平面上运动时[28]，通过对前后轮轮廓和水平面参数化处理，最后可以得到由广义坐标表示的两个几何约束方程

以上两个方程分别是广义坐标z和?关于广义坐标(θ,δ)的隐式函数，详细推导过程参见文献[28].显然，根据以上两个几何约束关系，z和?不再是独立的广义坐标，从而可以将自行车构型空间的维数由9维降为7维，且定义降维后的独立广义坐标向量为?q=[x,y,ψ,θ,δ,φr,φf]T.

假设自行车在水平面上作纯滚动，即每个接触点处速度为零.这一速度条件使得每个接触点处对应两个独立的速度约束：该接触点的速度在水平面上的两个分量为零.考虑自行车有两个接触点，故自行车系统共有4个纯滚动约束.

为恰当表示以上4个滚动约束方程，需要在接触点处建立局部的接触坐标系.该坐标系不一定要求是正交的，但要求局部接触坐标系中的一个坐标轴对应接触点处共切面的公法线方向[34].利用接触条件得到的接触体的轮廓参数关于广义坐标的函数关系，然后利用两个接触体在接触点处的相对速度的运动学关系，并将该相对速度为零的矢量条件映射到所建立的局部接触坐标系上.这样，在每个接触点处可以建立3个速度约束方程.容易证明，沿两个接触点处公法线方向的两个速度约束方程一定是两个几何约束方程(1)的全微分形式.自行车两个接触点引起的独立的速度约束方程仅有4个，它们对应4个独立的非完整约束方程[30].同时，这也意味着自行车系统的自由度为3，(这里的自由度采取分析力学中的传统定义，即系统的独立虚位移数[35]).

对受约束的多体系统来说，可以将几何约束和微分约束在速度空间中进行统一表示[36].对Whipple自行车系统来说，可以直接利用两个接触点相对速度为零的条件，并将以上条件映射到接触点的局部坐标系中，进而可以将相应的完整和非完整约束统一表示为

其中，是一个6×9的矩阵，其具体表达式见文献[28]的附录B.

当自行车在曲面上运动时，由于曲面接触点处的几何参数与描述自行车构型之间的广义坐标存在非常复杂的函数关系，通常需要采用计算机符 运算得到系统的约束方程.这一过程可以按照如下的流程进行:在两个接触点P和Q处分别定义三个曲面曲线参数，记为pi=[pi,1,pi,2,pi,3]T,i=r,f，其中，前两个曲面参数表征接触点在路面上的位置信息，最后一个曲线参数表征接触点在轮缘上的位置信息.根据文献[34]给出的接触几何条件，可以建立以上6个参数与系统广义坐标q之间的函数关系以及两个几何约束方程如下

结合以上关系，并利用接触点处相对速度等于零的条件，将接触点的相对速度表示在局部接触坐标系中，可以得到自行车在任意曲面上运动时的速度约束方程[28]

2、动力学建模

对非完整系统来说，采用Gibbs-Appell方法可以得到具有最小维数的动力学方程[19,32,45].Xiong等[28]基于约束方程选取了3个独立的准速度，然后利用Gibbs-Appell方法推导得到了自行车的约化动力学方程

其中，系数cii,k(q)是关于广义坐标q的函数，Pi(q)为广义力.通过数学软件MAPLE符 运算可以发现，方程(5)中的系数具有如下性质：

(1)各项系数只与σ1和σ2有关，即

当自行车在空间曲面上运动时，接触点处约束方程的数学描述更为复杂.这对动力学方程组的约化带来了很大的困难.Xiong等[27]利用第一类拉格朗日方程[49]得到了自行车在曲面路面的运动方程，对应为一组微分代数方程

其中，与自行车的重力和惯性力有关，是完整约束和非完整约束的系数矩阵，是前后车轮在接触点处与完整约束和非完整约束相对应的6个拉格朗日乘子.

3、自行车的相对平衡点及其稳定性

目前关于自行车的相对平衡点及其稳定性的研究主要有3种情况：

(1)自行车在水平面上的匀速直线运动；

(2)自行车在水平面上的匀速圆周运动；

(3)自行车在旋转曲面上的匀速圆周运动.

为研究以上动态运动的稳定性，首先需要通过系统的动力学方程，得到其对应的平衡解，然后在平衡位置附近建立线性化动力学模型，通过分析该线性系统的Lyapunov稳定性问题，确定系统相对平衡点的稳定性.历史上，人们对Whipple自行车做匀速直线运动情况的稳定性问题研究的较多[17,50,51]Carvallo[9]和Whipple[10]于1899年分别独立推导出自行车的线性化方程组，并针对具体的模型参数计算了自行车自稳定的速度范围.文献[11,15,17,42]相继针对Carvallo-Whipple自行车构型，分别采用Lagrange方法和矢量力学方法推导出了线性化方程.Schwab[42]与Meijaard[17]还利用多体系统动力学计算软件分别求出了线性化方程的系数，并与手算结果进行了对比验证.Lennartsson[52]和Astr¨om[53]利用多体动力学软件Sophia得到的非线性动力学方程在平衡位置附近进行线性化处理，也得到了相应的线性化方程组.Basu-Mandal[18]利用对带乘子的Lagrange方程线性化，Moore[25]对Kane方法得到的非线性动力学方程线性化.以上方法均得到了自行车做匀速直线运动时正确的线性化方程组，它具有如下的标准格式

其中，η=[θ,δ]T,f=[P1,P2]T，ω0是自行车匀速直线运动时后轮的角速度.其他系数矩阵M,K0,C和D是与自行车几何和惯量参数相关的常数矩阵.具体表达式见文献[28]的附录.

近些年来，学者们发现自行车的自稳定性与自行车的所有参数相关[26,54]，在经典Whipple自行车模型基础上还应进一步考虑非理想因素，如可以考虑车轮的滚动摩阻效应[55]；前车架系统对稳定性的影响[36,49,51,56,57,58]，其中前车架系统包括车头前立柱、前叉和前轮，在该模型中考虑了前车架系统的相对结构运动导致的水平速度自由度；自行车轮胎特征、路面梯度以及阻尼对稳定性也有影响[36,40,59,60,61]；自行车车轮的侧面滑动影响[11]；还有学者考虑了侧面风对自行车自稳定性影响[62].

自行车匀速直线运动的稳定性可以根据线性化动力学方程的特征根进行判断，Meijaard等[17]得到了Whipple自行车在运动速度0～10m/s范围内动力学方程(7)的特征根分布情况.根据特征根的分布，可以将自行车匀速直线运动受到小扰动时的状态分为weave运动，capsize运动和caster运动.其中，自行车稳定直线运动的速度限定在4.3m/s<v<6m/s.Doria等[35]考虑自行车的前车架系统，轮胎特征等影响，得到了4种不同动力学模型下的特征根的实部和虚部与速度的关系，并提出了局部稳定指标和整体稳定指标.

除此之外，也有学者对自行车在水平面上的匀速圆周运动进行了研究.Psiaki[13]和Franke等[63]分别得到了一族匀速圆周运动解，而Lennartsson[52]和Astr¨om等[53]则得到了三族匀速圆周运动解.BasuMandal等[18]首次给出了Whipple自行车在水平路面上的全部匀速圆周运动解族，他们一共得到了4个解族，而其中的两个解族在某种意义上可以合并为一个解族.该结果得到了Wang等[33]和Xiong等[28]的验证.

为了得到自行车在水平面上的匀速直线运动和匀速圆周运动的统一描述，并进行更加严格的稳定性分析，Xiong等[28]基于Gibbs-Appell方程对这两种运动进行了研究.这两种运动在某种意义上都可以看作是自行车动力系统的平衡点.注意到σ3=φr是循环坐标，可以将式(5)改写成一阶常微分方程组的形式

其中.设x0是方程(8)的平衡点，即其满足f(x0)=0，则必有x0=(θ0,δ0,0,0,ω0)，其中θ0,δ0,ω0是3个常数.根据前文提到的方程(5)系数所具有的性质，这3个常数只需要满足两个代数方程.因此，在这两个代数方程独立的情况下，这样的平衡点不是孤立的，而是以单参数解族的形式存在的.进一步，平衡点可以分为两类.第一类平衡点为，其中ω0为任意常数，这样的平衡点代表自行车的匀速直线运动.第二类平衡点的θ0,δ0不同时为零，这样的平衡点代表自行车的匀速圆周运动.利用Whipple自行车的经典参数[17],Xiong等[28]得到了四族单参数圆周解，与Basu-Mandal[18]及Wang等[33]的结果完全一致.

对于一些特殊的动力学系统[64]，其平衡点的Lyapunov稳定性可以通过构造Lyapunov函数，即采用Lyapunov直接方法[65]进行判断.但是，对于自行车系统，其动力学方程组系数的复杂性，导致相应的Lyapunov函数不易被求解.因此，Xiong等[28]基于Lyapunov间接方法，即根据系统(8)在平衡点处Jacobi矩阵的特征根分布来判断平衡点的Lyapunov稳定性，并且首次给出了自行车匀速直线运动和匀速圆周运动的非线性稳定性分析.事实上，由于x0是非孤立的，可以证明该点的Jacobi矩阵必有一个零特征根，这意味着x0是非双曲的平衡点.因此，x0的非线性稳定性与线性稳定性并不直接等价[66].Xiong等[28]结合中心流形定理[67]，第一次严格地证明了如果平衡点x0处的Jacobi矩阵其余4个特征根均具有负实部，则x0是稳定的，但不是渐近稳定的.

对于自行车在旋转曲面上匀速圆周运动的稳定性，Xiong等[27]亦开展了相关研究.首先，他们给出了自行车在旋转曲面上圆周解的定义.该定义给定了广义坐标q，曲面曲线参数p以及Lagrange乘子A随时间变化的规律，其中一共有21个待定常数，由微分代数方程组(6)确定.为了消除时变变量，他们引入了一系列旋转变化，将圆周解等价转化为了一个新的微分代数方程组的平衡解，其中一共有21个待定常数.同时他们发现，这21个常数只需要满足20个代数方程.因此与在水平路面上运动的情形一样，这样的平衡点一般来说不是孤立的，而是以单参数解族的形式存在的.他们以旋转抛物面ζ=α(ξ2+η2)，α0为对象，在不同的曲面系数下都找到了4个单参数圆周解族.这些圆周解族的结构随着旋转曲面参数的变化而演化，其中包含了不同解族之间的相互碰撞、分裂等现象.

旋转曲面上匀速圆周运动的稳定性指的是相应平衡解的Lyapunov稳定性.在平衡点处将新的微分代数方程组线性化，并利用线性约束方程消去非独立的变量，可以得到一个7维的线性系统.同样，平衡点的非孤立性意味着该线性系统必有一个零特征根.结合中心流形定理可以证明，如果线性系统的其余6个特征根都具有负实部，则该平衡点是稳定的，但不是渐近稳定的.数值计算结果发现，随着旋转抛物面系数α增大，稳定圆周解的区间越来越小.当α增大到某个临界值时，所有的圆周解都是不稳定的.

自行车的自稳定性研究除了基于解析模型的定量分析外，还可以从物理角度得到更直观的认识[26,68].历史上，对自行车稳定性的观点主要有3种：一是以Appell[69]、Timoshenko和Young[8]为代表的离心力观点，即前轮朝车体倾斜方向转动使车体转为曲线运动，所产生的离心惯性力平衡了重力的倾覆作用；二是以Klein和Sommerfeld[70]为代表的陀螺效应观点(如图2所示)，即车轮旋转产生的陀螺力矩在自稳定现象中起到了重要作用；三是以Grammel[71]为代表的“脚轮效应”观点(如图3所示)，即前轮尾迹c是影响自行车稳定性的关键因素，当c>0时，前轮的法向支撑力产生对前叉转轴的力矩，推动车把朝倾斜方向转动，使离心力效应的稳定作用自动实现.Jones[72]安装了一个与前轮并列反向旋转的副车轮，以消除陀螺力矩作用，用实验否定了陀螺效应，验证了脚轮效应.刘延柱[73]通过简化的力学分析对这三种效应都作出了解释.在另一篇文献[74]中，他经过估算发现在受控情况下，车轮的陀螺效应仅占稳定性因素的3%左右.

图2前轮的陀螺效应

图3车轮脚轮效应

但是，Kooijman等[26]则否认了陀螺效应和脚轮效应是保证自行车实现稳定性的必要条件.他们设计了一款“TBS自行车”(如图4所示)，该“自行车”同样包括车身、前叉、前轮和后轮4个刚体，但是车身和前叉简化为各自带有集中质量的杆件.在该实验自行车中，一方面通过在前轮附近利用反向旋转的副轮彻底消除陀螺效应，另一方面前叉转轴延长线与地面的交点位于前轮触地点的后方，与Jones[72]的自行车正好相反，以此消除脚轮效应.通过实验发现该“自行车”仍然可以实现自稳定，因此，他们认为陀螺效应和脚轮效应不是自行车实现自稳定的必要条件，但也不能否认这两个条件对自稳定的影响效应.他们同时验证了自行车的其他参数，如前车架质量分布，各个几何参数的相互耦合作用等等，均对自行车自稳定性有重要影响.刘延柱[75]认为Kooijman等[26]的自行车并非无脚轮效应，他通过理论分析发现，由于其重心高度远大于普通自行车，前轮的侧向摩擦力的正面作用超过法向支撑力的负面作用，由此推动前叉朝脚轮效应正确方向转动.

图4TBS自行车试验模型[26]

为了明确自行车各个参数对稳定性的影响程度，Xiong等[28]对自行车匀速直线运动状态下各物理参数对稳定性的敏感性进行了理论分析.他们发现轴距w越大，自行车在水平面上作匀速直线运动时的自稳定区域越小；前轮尾迹c越大，自稳定区域越大，这与Jones[72]的观点基本一致；前叉转轴倾角λ影响相对较小，意味着该参数对自稳定具有很好的鲁棒性；前框架和后框架质心在惯性系中的水平位置对自稳定具有显著影响，这与Kooijman[26]的观点一致；而对于自行车的车轮参数而言，增大前车轮半径，减小后车轮半径则有助于提高自行车的自稳定性；与陀螺效应有关联的前轮惯量参数对自稳定起相反作用.另外，Doria[35]一方面基于方差的敏感度分析方法对自行车不同几何参数对稳定指数的影响进行了分析，另一方面利用相关分析方法，根据不同的相关系数判断几何参数对自行车系统整体稳定指标的影响规律.Tak等[76]基于自行车线性动力学模型，研究了不同设计参数对自行车偏转速度和车把摇摆速度的灵敏度，结果发现在自行车的25个参数中有7个参数对稳定性影响明显.

4、自行车实验研究及控制

为了验证相关理论分析的正确性，许多学者开展了关于自行车的动力学实验研究.Roland和Lynch[77]实验分析了自行车的横向动力学以及轮胎特征参数，将实验结果与他们所提出的理论模型[78,79]进行对比验证；Miyagishi[80]设计了自行车机器人，用于实现无人骑行自行车的稳定性运动.

为研究自行车的自稳定特性，Kooijman等[81]在自行车上安装了陀螺仪装置及速度传感器(如图5所示)，测量自行车运动过程中的车轮角速度、车身偏航速度、车把转角以及后轮转角等相关运动学量，并反演自行车线性化方程的特征根，验证匀速直线运动时的稳定性理论结果.Escalona等[82]设计了类似的自行车动力学测试系统.通过测量自行车运动过程中关键点的位移、速度和加速度量，验证自行车多体动力学系统模型和的正确性和有效性.

近年来，人们也开展了一系列关于自行车控制的研究.Lowell和McKell[83]综合考虑了前文所述的自行车的3个稳定因素，将驾车人对前叉的控制用简单的比例规律替代，建立数学模型计算，却得出不稳定的结果.刘延柱[84]对控制规律作了一些修改，即考虑前叉转轴的阻尼因素，且控制规律中除车体倾斜的角度信息以外，还计入了角速度信息，这样计算结果即转变为渐近稳定.有的学者通过在自行车上增加外置平衡装置，如飞轮[85,86]和平衡块[87]实现自行车的稳定持续前行.Tanaka[88]以车把转角为控制变量，采用PD控制器研究自行车的稳定行驶.Keo[89]提出了outputzeroing控制器的方法实现自行车在速度为零时的直立状态的稳定性；Owczarkowski[90]和Anjumol[91]采用LQR控制方法实现自行车的直线跟踪控制；Toshiyuki和Zhang等[92]采用带扰动观测器的滑模控制方法研究运动稳定性，Chen等[93]采用滑模控制器研究自行车的车轮转动角跟踪控制问题.除以之外，有些学者[23,94,95,96,97,98]以电动摩托车为对象，开展了相关动力学的实验研究工作，为未来新型自行车的研制和实验提供了必要的借鉴资料.

图5自行车的实验测试系统[85]

5、结论与展望