假设您正在骑自行车，有人从旁边推了你一下。为了迅速让自己保持平衡，你会向同一方向转动自行车的把手防止跌倒。骑车的人出于本能会这样做，但神奇的是，自行车也能做到这一点。现代自行车设计中的自平衡功能，能在其运动失控时很好地保持平衡。让我们一起来看看如何在 COMSOL Multiphysics 中模拟这种效应。
自行车的自稳定性自安全自行车设计于 19 世纪 80 年代面世以来，现代自行车本质上并未发生太多变化。一个多世纪后的今天，科学家们仍在努力寻找使自行车保持自平衡的关键。换句话说，怎样才能让不受控制的自行车保持平衡？已有大量的研究论文通过解析方程解释了自行车的运动动力学。其中，最早的一篇代表性作品是由 Francis Whipple 撰写的一篇论文，他在该论文中推导出了专用于计算正在发生倾斜的、无人操控自行车动力学的一般非线性方程。
长期以来，人们普遍认为自行车能够保持稳定的因素有 2 个：前轮的陀螺效应和脚轮效应（或称为轨迹量），由前转向轴在前轮触地点前方与地面接触时产生。最近，来自代尔夫特理工大学和康奈尔大学的一个研究（请参见参考文献 3）发表了一篇关于Whipple 自行车模型运动的线性方程的全面综述。他们通过自己的研究成果展示了自行车的自稳定性。该研究表明，这种现象并非仅由某一因素引起，而是由陀螺效应和脚轮效应、自行车几何结构、速度及质量分布等多种因素共同作用，使不受控的自行车保持直立。
受这项研究启发，我们开发了一个多体动力学模型来演示无人操控自行车的平衡行驶情况。
自行车在不同时刻的运动状态。
自行车多体模型自行车教程模型演示了无人操控自行车在平坦表面上受到侧向力的扰动，在行驶方向发生倾斜的运动情况。我们对此进行了扩展分析，研究当前进速度和前转向轴倾角发生变化时自行车的自稳定性。
为了建立自行车模型，我们做出以下假设：
假设所有组件均为刚性。假设所有接头均无摩擦。建模时，假设自行车车轮与地面形成刀形触点。为车轮施加纯滚动约束。假设自行车在平坦表面移动。假设该自行车模型没有后座车架，通过在后车架上添加质量来定义。这辆自行车由4个刚性组件组成：后轮；包括无后座车架在内的后车架；包括车把在内的前车架；前轮。虽然车轮的厚度有限，但仍假设其与地面只有一个接触点。通过一个铰链关节使后车架与后轮相连接。并使用另一个铰链关节连接两个车架。该关节的轴构成自行车的转向轴。第三个铰链关节用于连接前轮和前车架。
使用 COMSOL Multiphysics 模拟的自行车几何结构。
模型中使用的自行车几何结构如上图所示。自行车各组件的惯性特性由质量、惯性矩和质心位置定义。重要的几何参数包括车轮半径、轴距 (wb) 及转向轴倾角 (st)。转向轴倾角可以控制轨迹量 (c)。
自行车的示意图。
如前所述，假设车轮与地面之间是纯滚动运动。在纯滚动的情况下，自行车触地点的速度为零。在该模型中，为了减少计算量，车轮滚动时不考虑任何接触，而是使用另一种多体动力学方法来模拟这种运动。为了描述车轮的偏航、倾斜及自转运动，我们创建了三个铰链及三个对应的无质量连杆。施加了三个无滑移约束来限制车轮在前进方向、面内横向及面外方向的滑移。由于无法对速度施加约束，所以我们将其转换为位移约束。
接下来，我们将描述设置这些约束的过程，相关内容在刚性滚轮教程模型中进行了详细说明。
为确保车轮纯滚动运动，需要三个约束来限制车轮在三个方向的滑移。
车轮模型，显示了三个约束方向的车轮模型。
这些约束为：
前进方向无滑移：
{frac{dold{u}}{dt}.old{e}_{2}=rfrac{dold{ heta}_s}{dt}}面内横向无滑移：
frac{dold{u}}{dt}.old{e}_{3}=rfrac{dold{ heta}_{l}}{dt}面外方向无滑移：
frac{dold{u}}{dt}.old{e}_{4}=0式中，old{e}_{2}、 old{e}_{3}和old{e}_{4}分别是瞬时前进方向（倾轴）、面内横向（旋转轴）及面外方向(old{e}_{4}=old{e}_{2} imesold{e}_{3})；frac{dold{u}}{dt}是重心的平移速度；r是车轮半径；frac{dold{ heta}_{s}}{dt}是旋转角速度；且frac{dold{ heta}_{l}}{dt}是倾斜角速度。
由于不能对速度施加约束，这些约束在时间上离散，通过以下方程实现：
(old{u}-old{u}_{p}).old{e}_{2}=r(old{ heta}_{s}-old{ heta}_{sp})(old{u}-old{u}_{p}).old{e}_{3}=r(old{ heta}_{l}-old{ heta}_{lp})(old{u}-old{u}_{p}).old{e}_{4}=0式中，old{u}_{p}、old{ heta}_{sp}和old{ heta}_{lp}分别是上一时刻的位移矢量、旋转角和倾角。
这些时间离散的无滑移约束需要之前时步的解。刚体的平移、旋转和上一时步中的瞬时轴均通过瞬态求解器中的全局方程和 上一个解节点来存储。
使用 COMSOL 模拟自平衡自行车的运动以一个转向倾角为 18° 的自行车为例进行分析。该自行车的初始前进速度为 4.6 m/s。1 s 后，通过在很短的时间内施加 500N 的力来扰动自行车。这个力使自行车在前进轴方向产生倾斜。
在第 1s 内，自行车以恒定速度沿初始前进方向移动。接着，侧向力使自行车开始倾斜。由于是无人驾驶自行车，因此不能人为干预使自行车保持平衡。那么，会发生什么呢？我们观察到，随着自行车开始倾斜，车把开始向同一侧转向。这种向倾倒方向的校正动作使自行车开始自我校直。
随着自行车继续沿不同的前进轨迹上前移，它开始朝反方向倾斜。由于转向运动稍微滞后，几乎与倾斜同时发生，因此倾斜幅度较小。这种前后倾斜和转向扰动持续进行，最终阻止了自行车倾倒。自行车直立向前移动，并获得略微增加的前进速度。倾角和转向角的振荡，以及角速度都会逐渐减弱并最终消失。
无人驾驶自行车在平坦面上的向前倾斜运动，箭头表示自行车的倾斜方向。
自行车的倾斜、转向旋转（左）及相对角速度（右）模拟结果。
稳定性分析到目前为止，我们已经知道自行车有自平衡能力。研究表明，影响自行车稳定性的因素不止一个。自行车的几何结构、质量分布及前进速度均对其稳定性产生影响。为理解这一点，我们进一步分析来研究两个参数的影响：初始前进速度和转向轴倾角。以上述自行车模型为参考构型，即转向轴倾角为 18°，向前行驶速度为 4.6m/s，并执行参数分析来研究这两个因素的影响。
初始前进速度变化在静止状态下，自行车无法保持直立。为了研究自行车速度对自稳定性的影响，我们以 1m/s 的步长在 2.6m/s 到 6.6m/s 范围内改变自行车的初始前进速度。在 2.6m/s 和 3.6m/s 时，自行车倾斜幅度非常大，不能保持稳定。在 5.6 m/s 时，自行车的倾斜速度减为零，但倾角是一个非零值。尽管这种构型能保持稳定，但自行车仍会略微倾斜，并做圆周运动。在 6.6m/s 时，倾角和转向角随时间的增加而变大，导致自行车运动不稳定。
不稳定稳定不稳定2.6 m/s3.6 m/s4.6 m/s5.6 m/s6.6 m/s在稳定情况下，初始前进速度为 5.6m/s（左）；在不稳定情况下，初始前进速度为 6.6m/s（右）。
转向轴倾角变化转向总成对于自行车的自稳定性至关重要。如果自行车不能转向（例如，车把被卡住），就无法抵抗倾斜而摔倒。对于这一点，控制轨迹量的转向轴倾角，对自行车的自稳定性起着重要的作用。
此次分析中，我们以 1° 的步长在 15° 到 21° 范围内改变转向倾角，从而改变轨迹量，来研究其对自行车稳定性的影响。当倾角为 15° 时，倾角和转向角继续随着时间的增加而变大，使这种构型变得不稳定。自行车在 16°到 19° 的范围内能够保持稳定，而超过该角度范围则不稳定。当转向倾角大于 19° 时，倾角和转向角会振荡，且振荡次数会随时间的增加而增加，导致自行车变得不稳定。
不稳定稳定不稳定倾角15°16°17°18°19°20°21°轨迹0.066 m0.0706 m0.0753 m0.08 m0.0848 m0.0896 m0.0945 m2种自行车个不稳定情况，转向轴倾角分别为 15°（左）和 21°（右）。
结束语在这篇博客中，我们描述了如何在 COMSOL Multiphysics 中使用多体动力学模块模拟无人操控的自平衡自行车的运动情况。文中演示了如何通过方程在刚性车轮上施加无滑移的约束，并将这些约束与自行车多体模型相耦合。然后，分析了初始前进速度和转向轴倾角对自行车自稳定性的影响。通过计算这些参数，我们发现，在一种构型下保持稳定的自行车在其他构型中会变得不稳定。
自行车自稳定性同时受多种因素的影响。文中的分析揭示了自行车的自平衡性与其朝倾斜方向转向的能力相关，得出的结论与之前的研究完全一致。
扩展阅读教程模型：刚性车轮的滚动教程模型：模拟自行车在平面上的运动Meijaard, Jaap P., Jim M. Papadopoulos, Andy Ruina, and Arend L. Schwab. “Linearized dynamics equations for the balance and steer of a bicycle: a benchmark and review.” In proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 463, no. 2084, pp. 1955-1982. The Royal Society, 2007.博客文章中: 在瞬态模拟中使用之前的解操作