基于MATLAB的珠子动力学分析

张汉松,左佳玉,吕 东,高兆辉

(大连海洋大学 1.食品科学与工程学院；2.信息工程学院；3.海洋科技与环境学院,辽宁 大连 116023)

牛顿第二定律是在惯性参考系下研究物体运动的变化情况，所谓惯性参考系是指相对于惯性系静止或作匀速直线运动的参考系。而非惯性参考系是指相对于惯性系做变速直线运动或曲线运动的参考系[1]。在非惯性系中引入惯性力后，牛顿运动定律的分析方式便可适用于非惯性系动力学分析，就可用来探讨非惯性系中物体的运动规律[2,3]。当珠子处于旋转的圆环内，相对地面惯性系做加速运动，珠子的运动情况属于非惯性系动力学分析。由于旋转体系中珠子运动的复杂性，牛顿第二定律已不再适用，目前还不存在能够广泛适用的运动状态分析模型，因而很难讨论影响珠子运动的参数。MATLAB 软件具有强大的计算和可视化功能，可以较为精确地模拟实验的结果，可避开复杂的理论推导以及条件苛刻的实验测量，使复杂、抽象的现象变得具体、直观。目前，利用MATLAB 软件对大学物理实验进行模拟仿真主要集中在光学实验[4-8]、声学实验[9]以及电磁学实验[10]等，对力学以及物体运动实验的模拟还比较少[11]。本文通过MATLAB仿真软件，首先模拟了理想状态与有摩擦力状态下，珠子角速度变化，得到了珠子临界角速度的条件，然后讨论了圆环半径、圆环转速、摩擦系数对珠子运动情况的影响。

1 基本原理及理论分析

1.1 科氏力

科里奥利力有些地方也称作哥里奥利力，简称为科氏力[12]，是对旋转体系中进行直线运动的质点由于惯性相对于旋转体系产生的直线运动的偏移的一种描述。科里奥利力来自于物体运动所具有的惯性。科里奥利力的计算公式如下：

(1)

1.2 欧拉法

欧拉法是常微分方程的数值解法的一种，其基本思想是迭代[13]。其中分为前进的EULER法、后退的EULER法、改进的EULER法。所谓迭代，就是逐次替代，最后求出所要求的解，并达到一定的精度，误差可以很容易地计算出来。

(2)

式(2)中的未知量yn+1可用欧拉公式预测计算，公式如下：

(3)

这就是改进的欧拉法。关于改进欧拉法的精度计算，利用半步长的泰勒级数展开式，由y″(xn+1)=y″(xn+h)=y″(xn)+hy‴(η)

(4)

比较上述式(3)和式(4)后，利用yn=y(xn)得到改进欧拉法的精度计算公式：

(5)

1.3 拉格朗日方程

对于完整系统用广义坐标表示的动力方程，通常是指第二类拉格朗日方程，是法国数学家J.-L.拉格朗日首先导出的。通常可写成：

(6)

式(6)中T为系统用各广义坐标qj和各广义速度q′j所表示的动能；qj为对应于qj的广义力；N(N=3n-k)为这完整系统的自由度；n为系统的质点数；k为完整约束方程个数。

1.4 动力学分析

图1为珠子的运动示意图。圆环绕着垂直于直径的轴旋转，珠子处于旋转的圆环内，让小珠子在环内凹槽中滚动。对珠子的运动进行动力学分析，分珠子和圆环之间无摩擦力和有摩擦力两种情况考虑，讨论珠子(用角度来描述)在环上凹槽中的运动情况。

图1 珠子的运动示意图

以圆环为参考系，珠子受到支持力、重力、摩擦力和离心力。在旋转环的内侧，为了便于比较，首先得到点状珠子的运动。这种处理包括非垂直旋转轴的情况，这会引起一些以前没有发现的共振现象。由拉格朗日方程得到以下方程组[14]：

(7)

式(7)中，T是动能，m是珠子的质量，R是圆环半径，r是珠子半径，ω是圆环转动的角速度，g是重力加速度，θ是珠子的位置，FN是珠子受到的支持力，K是珠子受到的科式力，f是珠子受到的摩擦力。

2 仿真研究与分析

2.1 无摩擦力珠子的运动情况

首先考虑理想情况，珠子和圆环之间无摩擦力，即摩擦系数μ为0。此时有如下方程：

(8)

式(8)中，ω是圆环转动的角速度，R是圆环半径，g是重力加速度，θ是珠子的位置。

用MATLAB进行虚拟仿真，用欧拉法进行迭代计算珠子的轨迹，给定角速度(ω=80，μ=0)得到珠子的运动轨迹图，如图2所示。

时间/s

2.2 有摩擦力珠子的运动情况

探索火星有限公司的首席执行官克里斯·卡伯里（Chris Carberry）正在努力琢磨，如何能让“疯子”和“火星人”一起通力合作，因为单单哪一派都不会有任何前途。也许我们能成功地选择去往月球——再去往火星。“我们正在寻找达成那一目标的方式，它不会将人类登陆火星一事耽搁上数十年之久。”他说道。除了研发协作，或许商业公司可以奔着登月目标而去，政府则聚焦于登陆火星的任务。

(9)

式(9)中，μ是摩擦系数，FN是珠子受到的支持力，K是珠子受到的科式力，f是珠子受到的摩擦力。

同样采用MATLAB模拟仿真，得到珠子的实际运动轨迹，如图3所示。

时间/s

由图3可知，珠子的运动为阻尼振动，最后在平衡位置附近停止。由于摩擦力的关系，珠子最终不会恰好停在平衡位置。如果忽略摩擦力，平衡位置比实际位置要高，当存在摩擦力时，虽然也会在某一位置停止运动，但是会低于无摩擦力情况下的位置高度。

2.3 无摩擦力时，角速度ω对珠子位置的影响

因为珠子在圆环中为滚动摩擦，所以摩擦系数可以设置的很小μ=0.1或者暂时忽略不计。当摩擦系数μ为0时，通过改变圆环不同转速，起初设定角速度ω很小，达不到珠子的转动临界速度，珠子一直在圆环底部。当圆环转速大于等于临界角速度时珠子会沿圆环滚动，不同的角速度珠子达到的最大高度不同，上升到最大高度的时间也不同。研究不同角速度ω下珠子运行一个周期(小球从圆环最低点到最大高度再到圆环最低点)所需要的时间，即珠子的角度随时间的变化曲线。用MATLAB进行模拟仿真，结果如图4所示。

时间/s

由图4可知，随着角速度的增大珠子停滞的最大角度在随之变大，而且珠子完成一个周期所需要的时间在随之变短。当角速度达到一定值时，珠子的最大停滞位置就不会再发生变化或者说是变化极小。

由图5可知，摩擦系数μ和圆环半径R一定时，圆环转动角速度ω越大，珠子运动的最大角度位置越大。当角速度达到一定值时，珠子的最大停滞位置不再发生变化或者变化极小。

角速度/rpm

2.4 摩擦因数μ对珠子运动的影响

运用控制变量法，假设摩擦因数μ对珠子的转动有影响，而圆环材质的改变对摩擦因数的影响微乎其微，通过MATLAB进行模拟仿真，主要通过改变摩擦因数的大小来分析摩擦因数对珠子在圆环上运动的影响，仿真结果如图6所示。

摩擦系数μ

由图6可知，随着摩擦因数的改变，珠子转动所需的临界角速度ω0也会增大。这是由于摩擦系数的增大，珠子转动所需要的向心力更大，那么就需要更大的角速度。所以，随着摩擦系数的增大，临界角速度也在增大。

2.5 圆环半径R对珠子运动的影响

3 结 语

基于上述分析，采用MATLAB软件仿真对珠子的动力学进行了分析。结果表明，圆环旋转角速度决定珠子最终的平衡位置，珠子初始位置θ和摩擦系数μ决定珠子开始运动的临界角速度ω0，当圆环角速度ω大于临界角速度ω0时，珠子开始运动，最终停在平衡位置。珠子的运动过程类似阻尼振动，珠子到达的最大角度位置和圆环角速度ω、摩擦系数μ和圆环半径R有关，μ、R一定时，ω越大，珠子运动的最大角度位置越大；ω、R一定时，μ越大，珠子运动的最大角度位置越小；ω、μ一定时，R越大，珠子运动的最大角度位置越小。特征时间与角速度ω及摩擦系数μ有关，摩擦系数μ一定时，角速度ω越大，特征时间越小；角速度ω一定时，摩擦系数μ越大，特征时间越大。